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Chapter  1 


Executive  Summary 


The  aim  of  this  research  effort  was  to  extend  earlier  work  on  two-dimensional  airfoils  undergoing 
high-intensity  unsteady  motions  dominated  by  leading-edge  vortex  (LEV)  shedding  to  finite- wing 
geometries.  In  the  earlier  work  (supported  by  AFOSR  grant  FA9550- 10- 1-0120,  Pis:  Gopalarath- 
nam,  Edwards,  and  01,  PM:  Dr.  Douglas  Smith),  an  integrated  theoretical,  computational,  and 
experimental  research  effort  was  used,  in  which  experiments  and  higher-order  computations  were 
used  to  develop  a  low-order  method  for  unsteady  aerodynamic  analysis  of  airfoils  and  plates  under¬ 
going  large-amplitude,  high-rate  motions.  Owing  to  the  fact  that  such  flows  are  dominated  by  LEV 
shedding  and  flow  separation,  they  are  well  outside  the  validity  of  classical  theoretical  methods. 

The  major  contribution  from  the  earlier  research  effort  was  in  identifying  the  importance  of  leading- 
edge  suction  in  governing  the  initiation,  formation,  and  termination  of  vortex  shedding  from 
rounded  leading  edges  of  unsteady  airfoils  in  two-dimensional  flow.  It  was  shown  that  the  value  of 
this  leading-edge  suction  at  any  time  instant  in  unsteady  flow  could  be  tracked  in  an  unsteady  thin 
airfoil  theory  using  an  inviscid  parameter  which  was  named  the  Leading-Edge  Suction  Parameter,  or 
LESP.  When  the  instantaneous  LESP  exceeds  a  critical  value,  LEV  shedding  occurs.  One  of  the  ma¬ 
jor  insights  from  that  research  effort  was  that  the  critical  LESP  is  dependent  only  on  the  airfoil  and 
Reynolds  number,  and  is  largely  independent  of  motion  kinematics  so  long  as  the  LEV  formation  is 
not  preceded  by  significant  trailing-edge  separation.  Typically,  LEV  formation  without  accompa¬ 
nying  trailing-edge  separation  occurs  in  high-rate  and  high-reduced-frequency  motions.  Thus,  for 
this  class  of  motions,  the  critical  LESP  for  a  given  airfoil  and  Reynolds  number  can  be  determined 
from  CFD  or  experiment  for  one  prototypical  motion,  and  can  then  be  used  for  any  other  high-rate 
motion,  including  arbitrary  combinations  of  pitch,  plunge,  and  surge.  This  insight  was  used  to 
augment  an  inviscid  unsteady  thin  airfoil  theory1  with  the  addition  of  a  discrete- vortex  method  to 
handle  intermittent  LEV  shedding.  The  resulting  method,  named  the  LESP-modulated  discrete 
vortex  method,  or  LDVM,  is  described  in  detail  in  Ramesh  et  al. 2  The  results  from  LDVM  show 
remarkably  good  agreement  in  forces  and  flow  fields  with  computational  fluid  dynamics  (CFD)  and 
experiments.  The  advantage  of  the  rapid  computational  capability  offered  by  a  low-order  method 
like  the  LDVM  was  illustrated  in  the  prediction  of  non-linear  aeroelastic  behavior  of  airfoils  having 
LEV  shedding.  As  described  in  Ramesh  et  al.,3  such  aeroelastic  analysis  requires  computation 
for  at  least  1,000  convective  times  to  establish  the  long-time  aeroelastic  behavior,  which  would  be 
prohibitive  with  any  high-order  CFD-like  method. 

Motivated  by  the  success  of  the  two-dimensional-flow  work,  the  research  in  the  current  effort  was 
aimed  at  extending  a  three-dimensional  inviscid  analysis  method  to  handle  vortex  shedding  from 
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rounded  leading  edges  of  finite  wings.  Systematic  studies  of  finite  wings  with  LEV  shedding  were 
analyzed  using  an  incompressible  Reynolds- Averaged  Navier  Stokes  (RANS)  CFD  method,  the 
results  of  which  were  used  for  the  development  of  the  low-order  method.  At  the  foundation  of 
the  low-order  finite- wing  analysis  method  is  a  traditional  unsteady  vortex  lattice  method  (UVLM) . 
In  the  first  part  of  the  current  work,  the  objective  was  to  assess  the  applicability  of  the  LESP 
concept  to  prediction  of  the  time  instant  and  spanwise  location  of  LEV  initiation  on  a  finite  wing 
undergoing  unsteady  motion.  The  UVLM  was  used  to  calculate  the  spanwise  variation  of  LESP 
at  every  time  step.  LEV  initiation,  from  low-order  prediction,  is  assumed  to  occur  at  the  time 
instant  and  spanwise  location  when  the  local  value  of  LESP  equals  the  two-dimensional  value 
of  critical  LESP.  Comparison  of  the  low-order  prediction  of  LEV  initiation  for  a  large  number 
of  wing  planforms  with  CFD  predictions  is  excellent,  confirming  that  the  flow  physics  governing 
LEV  initiation  in  finite  wings  is  the  same  as  that  for  airfoils.  Further,  the  critical  value  of  LESP 
obtained  for  two-dimensional  airfoil  flow  for  one  “high-rate”  motion  can  be  used  for  prediction  of 
LEV  initiation  on  any  finite- wing  geometry  undergoing  any  other  “high-rate”  motion.  Discrepancies 
between  the  low-  and  high-order  prediction  of  LEV  initiation  were  seen  in  very  low-aspect-ratio 
wing  (AR  =  2),  for  which  the  rolled- up  shed  vorticity  from  the  wing  tip  has  a  significant  influence 
on  the  wing  flow.  Since  the  UVLM  does  not  model  the  tip-vortex  roll  up,  there  is  a  noticeable 
discrepancy  for  the  AR  =  2  geometry. 

The  second  part  of  the  research  focused  on  modifying  the  UVLM  to  handle  LEV  formation.  In  this 
modified  UVLM,  a  vortex  sheet  from  the  leading  edge  is  modeled,  and  its  geometry  is  calculated 
in  each  time  step  by  assuming  that  the  corner  points  of  the  panels  forming  the  sheet  convect  with 
the  local  velocity.  The  resulting  self-induced  roll  up  of  the  sheet  is  simulated,  with  the  geometry 
agreeing  reasonably  well  with  CFD  results.  While  the  low-order  prediction  of  LEV  formation  shows 
promise,  significant  challenges  remain  in  the  modeling  of  the  sheet  roll  up  and  its  intersection  with 
the  wing  geometry.  Discrete-vortex  amalgamation  techniques,  studied  as  part  of  the  current  effort 
for  2D  flows,  shows  some  promise  for  reducing  the  number  of  vortex  lattices  in  LEV  sheets  from 
finite  wings.  It  is  likely  that  some  adaption  of  vortex  amalgamation  applied  to  vortex  sheets  from 
leading  edges  will  allow  for  robust  modeling  for  LEV  shedding  on  finite  wings. 

Suggested  future  work  includes  (i)  augmentation  of  the  current  UVLM  formulation  to  include  better 
modeling  of  tip-vortex  shedding,  (ii)  development  of  a  hybrid  vortex-sheet  /  vortex- filament  model 
for  LEV  shedding  from  finite  wings,  and  (iii)  extensions  to  the  UVLM  to  handle  LEV  formation, 
growth,  and  detachment,  allowing  for  modeling  of  intermittent  LEV  shedding. 
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Chapter  2 


Introduction 


Unsteady  aerodynamic  phenomena  are  prevalent  in  a  large  number  of  problems  in  modern  aerospace 
engineering  research.  These  include  dynamic  stall  in  wind  turbines  and  helicopter  rotors,  and 
flapping-wing  vehicle  (micro-air  vehicle)  design.  These  problems  are  characterised  by  apparent- 
mass  effects,  flow  separation  and  leading-edge  vortex  (LEV)  formation  and  shedding.  Unsteady 
flows  with  intermittent  LEV  formation  are  the  focus  of  the  current  research.  The  LEV  influences 
the  flowfield  tremendously,  and  is  responsible  for  an  enhancement  in  lift  while  it  is  present  above 
the  wing  and  also  large  nose-down  pitching  moments  and  flow  separation  over  the  entire  airfoil 
when  it  convects  either  off  the  trailing  edge  or  away  from  the  airfoil.1,4  In  examining  past  work 
to  provide  context  to  the  current  research,  we  first  examine  in  Section  2.1  LEV-dominated  two- 
dimensional  airfoil  flows  first  with  special  attention  to  low-order  modeling  using  discrete-vortex 
approaches.  Subsequently,  in  Section  2.2,  we  examine  relevant  literature  for  LEV-dominated  finite- 
wing  unsteady  aerodynamics. 


2.1  LEV-dominated  airfoil  flows 

Methods  of  simulating  the  physics  and  effects  of  unsteady  aerodynamic  phenomena  date  back  to 
Wagner5  and  Theodorsen.6  The  application  of  these  theoretical  methods  is  however  limited  by 
their  assumptions  such  as  attached  flow,  small  amplitude  motion  and  planar  wake  structure.  For 
unsteady  flows  with  vortex  shedding  such  as  those  considered  in  the  present  study,  experiments 
and  high-fidelity  computations  have  facilitated  fundamental  studies  of  the  underlying  phenomena. 
McGowan  et  al.,7  01  et  ah, 8  Garmann  V  Visbal9  and  Granlund  et  al.10  have  analysed  the  forces 
and  flowfields  for  such  unsteady  motions  over  a  broad  parameter  space  using  both  experimental 
and  computational  methods.  Pitt  Ford  &  Babinsky,11  Baik  et  al.12  and  Rival  et  al.13  have  studied 
leading  edge  vortices  using  experimental  techniques.  These  methods  are  not  suitable  for  the  initial 
phases  of  aerodynamic/control  design  because  of  cost  and  time  considerations.  This  has  been  the 
motivation  of  many  researchers  to  construct  low-order  models  for  unsteady  simulation  of  wings  and 
airfoils. 

Discrete-vortex  methods  have  been  successfully  used  in  past  decades  to  model  unsteady  separated 
flows.  These  methods  are  usually  based  on  potential- flow  theory,  and  the  shear  layers  depicting 
separated  flow  are  shed  from  the  surface  in  the  form  of  discrete  vortices.  Clements  &  Maull14 
and  Saffman  V  Baker15  have  written  detailed  reviews  on  the  historical  development  of  the  discrete- 
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vortex  method.  A  review  on  more  recent  applications  of  vortex  methods  for  flow  simulation  is  given 
by  Leonard.16  Sarpkaya,17  Clements,18  Kiya  &  Arie,19  amongst  other  researchers,  have  applied 
this  category  of  methods  successfully  to  model  flow  past  inclined  plates  and  bluff  bodies.  Katz20 
has  developed  a  discrete-vortex  method  for  separated  flow  past  an  airfoil,  where  the  location  of 
separation  on  the  airfoil  has  to  be  prescribed  using  information  from  experiments  or  other  means. 
More  recently,  low-order  methods  based  on  discrete  vortices  have  been  developed  by  Ansari  et 
ah, 21  Wang  &  Eldredge,22  Ramesh  et  al.2  to  model  leading  edge  vortices  in  unsteady  flows,  with 
applications  toward  insect  flight  and  MAV  aerodynamics.  Though  these  methods  are  based  on 
potential  theory,  they  capture  the  essential  physics  in  flows  of  interest  by  combining  inviscid  theory 
with  discrete-vortex  shedding.  Apart  from  providing  a  means  to  calculate  the  force  coefficients  on 
the  airfoil,  these  methods  also  enable  the  study  of  flow  features.  These  are  significant  advantages 
of  this  class  of  methods  over  semi-empirical  methods,  which  only  allow  determination  of  the  force 
coefficients  through  empirical  fitting. 

Many  of  the  methods  cited  above  assume  some  ad-hoc  start  and  stop  criteria  for  discrete-vortex 
shedding,  such  as  continuous  shedding  from  a  given  location  (valid  only  for  sharp  edges)  or  shedding 
that  starts  and  stops  depending  on  whether  the  local  angle  of  attack  exceeds  a  critical  value 
(valid  only  for  a  small  range  of  motions).  A  more  general  vortex  shedding  criterion  is  required  to 
make  discrete- vortex  methods  broadly  applicable  to  a  wide  range  of  geometries  (including  airfoils 
with  rounded  leading  edges)  and  arbitrary  unsteady  motions.  Ramesh  et  al.2  have  developed  a 
discrete- vortex  aerodynamic  method  to  model  unsteady  flows  with  intermittent  LEV  shedding  using 
a  leading-edge  suction  parameter  (LESP).  The  unique  aspect  of  this  method  (LESP-modulated 
discrete  vortex  method,  or,  LDVM)  is  that  vortex  shedding  is  turned  on  or  off  at  the  leading  edge 
using  a  criticality  condition.  This  method  is,  therefore,  ideally  suited  to  modeling  oscillatory  airfoil 
flows  in  which  intermittent  LEV  shedding  is  a  key  feature.  In  comparison  with  semi-empirical 
methods  in  which  several  parameters  are  typically  used,  for  a  given  airfoil  and  a  given  Reynolds 
number,  this  model  uses  only  a  single  empirical  constant,  the  critical  LESP,  and  is  a  highly  physics- 
based  approach. 

LDVM  has  been  shown  to  be  successful  in  predicting  forces  and  moments  on  the  airfoil,  as  well 
as  flow  field  around  the  airfoil  for  high-frequency  unsteady  maneuvers.  However,  it  is  necessary 
to  track  a  significant  number  of  discrete  vortices  in  order  to  obtain  the  instantaneous  forces  and 
moments  on  the  airfoil.  The  computational  complexity  increases  as  0(n2)  (when  fast  summation 
methods  are  not  used),  where  n  is  the  number  of  vortices  in  the  flowfield,  resulting  in  possibly 
large  computing  times.  The  computational  time  of  any  discrete  vortex  model  can  be  significantly 
decreased  by  reducing  the  number  of  discrete  vortices.  The  current  study  focuses  on  obtaining  a 
model  with  a  reduced  number  of  leading  edge  vortices,  thus  improving  the  computation  time. 

The  popular  Brown-Michael  model23  was  one  of  the  earliest  efforts  towards  modeling  the  leading 
edge  vorticity  using  reduced  number  of  vortices.  Specifically,  the  Brown-Michael  model  used  a 
single  vortex  with  time-varying  strength  to  represent  the  vorticity  shed  from  a  delta  wing.  Wang 
and  Eldredge22  recently  improved  on  this  model  to  derive  analytic  expressions  for  the  evolution  of 
single  vortices  using  impulse  matching  for  a  flat  plate.  Howe24  has  presented  a  generalized  correction 
for  the  Brown-Michael  model,  and  uses  the  model  to  study  the  effect  of  the  translating  vortex  of 
time-varying  strength  on  a  rigid  half-plane.  One  vortex  with  time-varying  strength  models  the 
leading  edge  vorticity,  while  another  similar  vortex  is  used  to  represent  the  trailing  edge  vorticity. 
Cortelezzi  and  Leonard25  used  a  single  vortex  with  time-dependent  strength  to  model  the  shear- 
layer  roll  up  of  a  semi-infinite  plate.  This  model  was  further  revised  to  accommodate  a  shear  layer 
feeding  vorticity  into  a  growing  vortex.26  Some  of  these  models  use  an  unsteady  Kutta  condition 
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at  the  leading  edge  to  find  the  strength  of  the  single  vortex  that  models  the  LEV. 

Several  experimental  studies  have  been  conducted  by  various  researchers  to  understand  the  evolu¬ 
tion  of  the  leading-edge-vortex  structure.  The  shear  layer  emanating  from  the  leading  edge  starts 
rolling  up  into  a  concentrated  vortex.  Vorticity  is  fed  into  this  concentrated  vortex  structure  by 
the  shear  layer,  and  it  grows  in  strength.  In  case  of  prolonged  vorticity  shedding,  the  leading  edge 
vortex  pinches  off  from  the  shear  layer,  and  is  convected  downstream.  Meanwhile,  a  new  vortex 
roll-up  is  initiated  near  the  leading  edge.13,27,28  Antonini  et  al.29  use  semi-empirical  data  to  model 
the  expanse  of  the  time- varying  vortex  core. 

Several  articles  can  be  found  in  the  literature  that  address  the  merging  of  discrete  vortices  to  reduce 
the  count  of  discrete  vortices  in  the  flow  field. 17,30,31  For  example,  Nair  and  Taira32  give  a  network 
theoretic  approach  to  identify  important  vortex- vortex  interactions,  and  obtain  a  sparsified  model 
that  accurately  predicts  the  dynamics  of  the  original  system  based  on  these  interactions.  The 
conditions  used  for  merging  are  based  on  closeness,  relative  velocity  etc.  between  vortex  pairs. 


2.2  LEV-dominated  finite- wing  flows 

Vortex  shedding  from  the  leading  edges  of  wings  and  rotor  blades  have  been  observed  in  nature — 
on  swimming  and  flying  animals  and  seeds,33-36  and  in  engineering — on  rotorcraft,37  wind  tur¬ 
bines,38  swept  and  delta  wings,39,40  and  micro-air  vehicles41  and  flapping- wing  energy- harvesting 
devices.42,43  Extensive  investigation  of  leading-edge  vortex  (LEV)  formation  and  shedding  from 
airfoils  in  two-dimensional  flow  have  revealed  the  connection  between  the  onset  of  LEV  formation 
and  flow  separation  at  the  leading  edge  or  leading-edge  stall,  and  their  the  dependence  on  leading- 
edge  radius,37,44  Reynolds  number,37,44,45  and  unsteady  motion  kinematics  of  the  airfoil.10,46,47 
Numerous  computational  and  experimental  studies  have  shown  the  effects  of  the  LEV  growth,  po¬ 
sition,  and  detachment  on  the  forces  and  moments  experienced  by  the  airfoil.47  49  The  restriction 
to  two-dimensional  flow  also  enables  the  use  of  low-order  discrete-vortex  methods  for  modelling 
the  airfoil  LEV  formation  and  its  effects.50,51  In  contrast,  LEV  formation,  shedding,  growth,  and 
their  effects  on  finite  wings  are  considerably  more  complicated  by  the  presence  of  spanwise  ve¬ 
locities  and  pressure  gradients,  interaction  with  root  and  tip  vortex  structures,  and  the  interplay 
between  vorticity  production  and  spanwise/chordwise  advection. 52,53  On  some  geometries  like  the 
highly-swept  leading  edges  of  delta  wings,  vorticity  production  at  the  leading  edges  is  balanced  by 
spanwise  vorticity  transport,  leading  to  body-relative  stable  or  stationary  vortex  structures,  which 
can  be  harnessed  for  lift  enhancement  at  high  angles  of  attack  and  extra  maneuverability.54  The 
effects  of  these  LEV  flows  are  also  amenable  to  simple  and  elegant  theories  such  as  the  Polhamus 
leading-edge  suction  analogy.55,56 

On  other  configurations  like  unswept  wings,  the  absence  of  mechanisms  for  spanwise  transport 
of  shed  leading-edge  vorticity  appears  to  be  the  cause  for  non-uniform  shedding  and  chordwise 
advection,  leading  to  interesting  and  important  flow  structures  like  the  omega-shaped  (or  horse- 
shoe-shaped)  vortical  structures  that  have  been  observed  in  experiments  and  computations.57-59 
Further  complicating  the  three-dimensional  situation  are  the  rotational  effects  on  these  phenomena 
on  rotor  blades  and  flapping  wings.60  It  may  be  argued  that  the  first  step  in  unravelling  the  flow 
physics  of  finite- wing  LEVs  is  the  initiation  of  LEV  formation:  for  any  given  wing  and  motion 
kinematic,  at  what  time  or  angle  of  attack  (AoA)  and  where  along  the  span  does  the  LEV  start 
forming?  The  initial  portion  of  this  research  was  focussed  on  answering  this  specific  question. 
Building  on  an  earlier  work  on  initiation  of  LEV  formation  on  rounded-leading-edge  airfoils,2  and 
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using  results  from  three-dimensional  CFD  computations  for  a  large  number  of  finite  wings,  it  is 
shown  that  criticality  of  leading-edge  suction,  which  governs  LEV  formation  on  airfoils,  can  also  be 
reliably  used  to  predict  the  initiation  of  LEV  formation  on  finite  wings  in  low-Mach  number  flows. 


2.3  Layout  of  this  report 

Computational  (CFD)  and  experimental  studies  have  been  used  to  guide  the  development  of  the 
low-order  methods  developed  in  the  current  research.  Chapter  3  presents  the  CFD  methods  and 
Chapter  4  provides  a  description  of  the  experimental  facilty  and  techniques.  In  Chapter  5,  we 
present  a  detailed  study  of  LEV  initiation  in  airfoil  flows,  with  emphasis  on  determining  the  range 
of  validity  of  the  LESP  concept.  It  is  shown  that,  so  long  as  LEV  formation  is  not  preceded  by 
significant  trailing-edge  flow  reversal,  the  LESP  concept  holds.  This  condition  is  satisfied  for  high- 
rate  motions.  Chapter  6  presents  research  into  amalgamation  of  discrete  vortices  in  2D  airfoil  flows 
in  an  attempt  to  reduce  computational  time.  Ideas  learned  from  this  study  were  subsequently  used 
in  amalgamation  of  vortex  lattices  in  finite- wing  LEV  sheets  (Chapter  9).  The  main  focus  of  the 
current  effort  is  on  low-order  prediction  of  LEV-dominated  finite- wing  flows.  Chapter  7  presents  the 
unsteady  vortex  lattice  method  (UVLM),  which  forms  the  foundation  for  the  finite-wing  research. 
Also  discussed  are  the  modifications  to  the  UVLM  for  modeling  LEV  vortex  sheets. 

The  criterion  to  identify  rollup  in  the  model  presented  in  the  current  work  is  based  on  angular 
velocity  between  vortex  pairs  at  the  tip  of  the  shear  layer.  Additionally,  a  discrete  vortex  is 
merged  with  the  growing  vortex  based  on  its  closeness,  relative  and  angular  velocities  with  respect 
to  the  growing  LEV.  Chapter  8  presents  the  results  on  LEV  initiation  on  finite  wings,  including 
comparison  between  low-order  and  CFD  predictions.  Chapter  9  presents  results  from  the  modified 
UVLM  on  prediction  of  LEV  formation  on  finite  wings.  The  report  then  wraps  up  with  conclusions 
and  suggestions  for  future  work  in  Chapter  10. 
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Chapter  3 


Computational  Fluid  Dynamics 
Methods 


In  this  effort,  the  computational  fluid  dynamics  (CFD)  portion  of  the  research  had  two  objec¬ 
tives.  The  first  objective  was  to  provide  high-order  results  from  CFD  simulations  using  “state- 
of-the-practise”  Reynolds- Averaged  Navier-Stokes  (RANS)  simulations  and  Spalart-Allmaras  (SA) 
turbulence  model  to  guide  the  development  of  the  low-order  method  and  proide  validation  data. 
Section  3.1  provides  a  brief  description  of  this  portion  of  the  CFD  effort.  The  second  objective  of 
the  CFD  effort  was  to  develop  improved  capability  via  the  development  of  a  hybrid  large-eddy  sim¬ 
ulation  /  Reynolds- averaged  Navier-Stokes  method  suitable  for  3D  simulations  on  moving  meshes. 
Progress  on  this  portion  of  the  CFD  effort  is  described  in  Section  3.4.  A  version  of  this  method 
includes  the  Menter-Langtry  transition  model.  This  new  capability  is  expected  to  be  useful  for 
exploring  the  effects  of  transitional  boundary  layers  on  the  formation  of  LEVs  and  improved  pre¬ 
dictions  of  LEV-dominated  flows. 


Figure  3.1:  Body- fitted  mesh  used  in  computations. 
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Figure  3.2:  Vorticity  scale  for  all  CFD  vorticity  plots  in  this  report. 


3.1  RANS  CFD  for  supporting  the  development  of  the  low-order 
methods 

CFD  calculations  using  a  body  fitted  grid  were  performed  using  NCSU’s  REACTMB-INS  code, 
which  solves  the  time-dependent  incompressible  Navier-Stokes  equations  using  a  finite-volume 
method.  The  governing  equations  are  written  in  arbitrary  Lagrangian  /  Eulerian  (ALE)  form, 
which  enables  simulations  of  the  flow’s  response  to  the  motion  of  a  body-fitted  computational  mesh 
in  accord  with  prescribed  rate  laws.  Spatial  discretization  of  the  inviscid  fluxes  uses  a  low-diffusion 
flux-splitting  method  valid  in  the  incompressible  limit.61  This  method  is  extended  to  higher-order 
spatial  accuracy  using  TVD  interpolations  of  the  primitive  variables  [p,  u,  u,  re,  z>]T.  Viscous  terms 
are  discretized  using  second-order  central  differences.  A  dual-time  stepping  method  is  used  to  in¬ 
tegrate  the  equations  in  time.  An  artificial  compressibility  technique,  discretized  in  a  fully  implicit 
fashion  and  solved  approximately  using  ILU  decomposition,  is  used  to  advance  the  solution  in 
pseudo-time.  Typically,  eight  sub-iterations  per  physical  time  step  were  need  to  reduce  the  residual 
errors  two  orders  of  magnitude.  The  Spalart-Allmaras  model62  as  implemented  by  Edwards  and 
Chandra,63  is  used  for  turbulence  closure.  The  computations  were  performed  on  a  2-D  body-fitted 
mesh  containing  92400  cells  (Figure  3.1). 


3.2  Identification  of  LEV  formation  from  CFD 

The  surface  skin  friction  (Cf)  distributions  from  CFD  provide  a  quantitative  method  of  identifying 
the  onset  of  LEV  formation.  Formation  of  LEV  is  preceded  by  the  formation  of  a  small  region  of 
reversed  flow  with  negative,  or  counter-clockwise,  vorticity  at  the  surface  near  the  leading  edge  of 
the  airfoil.  Onset  of  LEV  shedding  is  subsequently  initiated  by  the  formation  of  a  shear  layer  at  the 
leading  edge  and  this  is  accompanied  by  the  development  of  a  small  region  of  positive,  or  clockwise, 
surface  vorticity  within  the  region  of  negative  vorticity  described  earlier.  Thus  the  formation  of 
alternating  negative  and  positive  vorticity  near  the  leading  edge  is  a  useful  signature  which  we  will 
use  to  determine  the  exact  instant  when  LEV  formation  begins.  Since  the  process  of  identifying 
onset  of  LEV  formation  visually  using  the  vorticity  plots  is  qualitative,  we  translate  the  same 
criterion  to  the  skin  friction  coefficient.  The  formation  of  a  negative  vorticity  region  at  the  surface 
near  the  leading  edge  is  indicated  by  a  negative  skin-friction  coefficient  near  the  leading  edge.  The 
formation  of  the  positive  vorticity  region  at  the  surface,  which  accompanies  the  development  of 
the  shear  layer,  can  be  identified  by  alternating  positive  and  negative  spikes  in  the  skin-friction 
coefficient.  Using  the  vorticity  scale  shown  in  Fig.  3.2,  Fig.  3.3  shows  a  series  of  vorticity  and 
skin-friction  coefficient  plots  for  an  airfoil  experiencing  LEV  formation. 

In  Figure  3.3,  the  vorticity  and  the  skin- friction  coefficient  plots  at  the  leading  edge  (from  x/c  —  0 
to  x/c  =  2%)  are  shown  during  four  instants  for  an  airfoil  undergoing  LEV  formation.  The  flow 
behavior  at  instants  (a)  through  (d)  in  the  figure  is  described  as  follows: 
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Figure  3.3:  Vorticity  and  skin  friction  coefficient  plots  from  CFD  during  the  LEV  formation  process. 
Appearance  of  positive  vorticity  region  is  circled. 


•  (a)  The  flow  is  attached  at  the  leading  edge. 

•  (b)  The  boundary  layer  is  still  attached.  However,  there  is  a  region  of  reversed  flow  near  the 
leading  edge,  reflected  in  both  the  vorticity  and  the  skin-friction  coefficient  plots. 

•  (c)  Start  of  formation  of  the  shear  layer:  A  small  region  of  positive  vorticity  starts  developing 
at  the  surface  within  the  negative  vorticity  region  at  the  leading  edge.  This  behavior  also 
corresponds  to  the  first  occurrence  of  a  small  region  of  positive  Cf  within  the  region  of 
negative  Cf.  In  the  Cf  plot,  the  formation  of  spikes  reaching  up  to  zero  and  positive  values 
in  the  region  of  negative  Cf  near  the  leading  edge  is  seen.  We  use  this  signature  in  the  Cf 
to  quantitatively  identify  the  time  instant  corresponding  to  the  start  of  LEV  formation  in  all 
results  from  CFD.  This  identification  is  done  by  post-processing  the  CFD  results. 

•  (d)  Shear  layer  is  well  established  and  LEV  shedding  is  in  progress. 

3.3  RANS  analysis  of  LEV  formation  on  finite  wings  and  compar¬ 
ison  to  experiment 

NCSU’s  REACTMB-INS  solver  is  used  for  the  computational  fluid  dynamics  (CFD)  calculations 
performed  in  this  study.  This  finite-volume  solver  formulates  the  time- dependent  incompressible 
Navier-Stokes  equations  in  an  arbitrary  Lagrangian/Eulerian  (ALE)  fashion.  The  ALE  form  enables 
moving-mesh  flow  simulations  on  the  3-D  body-fitted  computational  mesh.  An  incompressible 
version  of  Edwards’  Low-Diffusion  Flux  Splitting  Scheme  (LDFSS)61  is  used  for  discretizing  inviscid 
fluxes  in  space.  Discretization  of  viscous  terms  is  performed  using  a  second-order  central  difference 
method.  The  LDFSS  method  is  extended  to  higher  order  of  accuracy  in  space  using  the  Piecewise- 
Parabolic  Method  (PPM).64  For  time  integration,  second-order  temporal  accuracy  is  achieved  by 
using  an  implicit  artificial  compressibility  method61  with  subiterations  at  each  physical  time  step 
for  continuity  equation  convergence.  A  version  of  Spalart-Allmaras  one  equation  eddy  viscosity 
model,  modified  by  Edwards  and  Chandra,63  is  used  for  turbulence  closure. 


9 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Figure  3.4:  Representative  mesh  distribution  for  CFD  analysis. 


Figure  3.4  shows  the  mesh  distribution  of  a  rectangular  half- wing  used  for  the  finite- wing  cal¬ 
culations  in  the  current  work.  The  chord  length  is  0.1  meters.  The  O-type  mesh  has  164  cells 
chord- wise,  with  finer  resolution  near  the  leading  edge  and  trailing  edge.  The  span- wise  average 
spacing  on  the  airfoil  is  chord/100,  with  finer  resolution  near  the  tip  of  the  wing.  The  span- wise 
calculation  domain  extends  2  chord  lengths  beyond  the  tip  of  the  wing,  with  an  average  spacing  of 
chord/40  in  this  region.  In  the  wall- normal  direction,  cell  spacing  starts  as  5  x  10-5  meters  next  to 
the  wall,  and  has  a  growth  factor  of  1.15  moving  away  from  the  surface  until  the  spacing  reaches 
chord/100.  From  there,  cell  spacing  is  kept  nearly  uniform  of  chord/100  up  to  1.3  chord  from  the 
surface.  Then  coarser  meshes  with  a  growth  factor  of  1.15  extend  to  12  chord  lengths  away  from 
the  wing  surface.  Only  the  rectangular  wing  with  an  aspect  ratio  of  6  is  shown  here,  but  the  general 
guidelines  above  are  applied  to  meshes  of  all  other  wing  geometries  considered  in  this  study. 

The  CFD  model  is  validated  by  comparing  the  flow  solution  with  the  PIV  results  from  the  ex¬ 
perimental  study  of  Yilmaz  and  Rockwell59  for  a  rectangular  flat  plate  with  an  aspect  ratio  of  2 
undergoing  a  0-to-45-degree  pitch-up  motion.  Figure  3.5  compares  predicted  iso-surfaces  of  the 
second- invarient  of  the  velocity  gradient  tensor  (Q=5)  with  experimental  images  obtained  from 
the  PIV  database.  Side  views  of  the  3D  streamline  patterns  at  five  instances  in  time  are  shown 
in  Figure  3.6.  Compared  with  experimental  data,  the  streamline  patterns  from  CFD  simulation 
show  the  same  stage  of  development  of  the  LEV  at  each  time  instance.  Overall,  the  CFD  results 
compare  well  with  the  PIV  results,  giving  confidence  in  the  utility  of  the  CFD  technique  for  the 
present  work. 


3.4  LES/RANS  Modeling  Activities 

3.4.1  Model  Formulation 

As  some  three-dimensional  calculations  were  planned  for  this  work,  efforts  were  directed  toward 
developing  a  hybrid  large-eddy  simulation  /  Reynolds- averaged  Navier-Stokes  method  suitable  for 
3D  simulations  on  moving  meshes.  In  the  present  study,  the  transition  between  a  RANS  component 
(used  very  near  solid  surfaces)  and  the  LES  component  (used  in  the  outer  parts  of  developing 
turbulent  boundary  layers  and  in  free  shear  layers)  is  facilitated  by  the  action  of  a  flow-dependent 
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Figure  3.5:  Volumes  of  iso-Q  as  a  function  of  angle  of  attack.  Left:  PIV  experiment.  Right:  CFD 
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Figure  3.6:  Side  views  of  three-dimensional  streamline  patterns  as  a  function  of  angle  of  attack. 
Left:  PIV  experiment.  Right:  CFD  simulation. 
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blending  function,  which  modifies  the  eddy  viscosity  field  as  follows: 


Pf  =  P 


(1  ^)^t,sgs 


UO 


(i) 


where  Y  is  a  time-dependent  blending  function  that  connects  the  RANS  and  LES  branches.  A 
model  due  to  Lenormand,  et  al.65  is  currently  used  for  the  subgrid-scale  eddy  viscosity.  The 
blending  function  is  generally  designed  to  transition  the  model  from  RANS  to  LES  approximately 
as  the  boundary  layer  shifts  from  its  logarithmic  to  its  wake-like  structure.  As  such,  the  RANS 
component  acts  as  a  wall-layer  model  for  the  majority  of  the  flow,  which  is  modeled  as  a  large-eddy 
simulation.  The  most  recent  model66,67  introduces  an  outer-layer  turbulence  length  scale  into  the 
argument  for  T  as  a  means  of  removing  a  calibration  procedure  required  in  an  earlier  version.68 
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In  this  expression,  k  is  the  ensemble-averaged  modeled  turbulence  kinetic  energy,  kf>  is  the  ensemble- 
averaged  resolved  turbulence  kinetic  energy,  uo  and  uo  are  instantaneous  and  ensemble-averaged 
modeled  turbulence  frequencies,  and  d  is  the  distance  to  the  nearest  wall.  The  combination  of 
instantaneous  and  ensemble- averaged  data  allows  the  RANS-to-LES  transition  position  Y  =  1/2  to 
fluctuate  about  a  mean  value  that  is  a  function  of  the  local,  ensemble-averaged  state  of  the  flow.  As 
it  is  dependent  on  both  inner-layer  and  outer-layer  turbulence  length  scale  information,  this  model 
is  more  capable  of  adjusting  to  departures  from  local  equilibrium,  and  a  problem- specific  selection 
of  a  model  constant  is  not  required.  The  required  ensemble  averages  are  currently  computed  using 
an  exponentially-weighted  moving  average. 

Also  tested  in  this  study  is  an  earlier  LES/RANS  model  developed  by  Choi,  et  al.68  This  model 
differs  from  the  above  in  that  a  pre-calibration  is  needed  to  calculate  a  model  constant  that  multi¬ 
plies  !  ^ —  in  Eq.  3.  This  model  constant  is  related  to  the  wall-coordinate  distance  in  which  an 

cfi  Kd^/u 

attached  boundary  layer  transitions  from  its  logarithmic  structure  to  its  wake-like  structure,  and 
it  can  be  computed  as  a  function  of  surface  distance  based  on  knowledge  of  the  boundary  layer 
thickness  and  edge  conditions.68  In  addition,  a  version  of  REACTMB  that  couples  the  LES/RANS 
formulation  of  Eq.  3  with  the  Menter-Langtry  transition  model69  was  developed.  This  model  re¬ 
quires  the  solution  of  two  additional  transport  equations,  one  associated  with  a  Reynolds  number 
based  on  momentum  thickness  and  the  other  associated  with  turbulent  intermittency.  As  the 
Menter-Langtry  model  only  alters  the  RANS  component  of  the  LES/RANS  model,  the  coupling 
was  straightforward  and  is  not  described  further  in  this  report. 
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3.4.2  Aerospatiale  A-Airfoil 


The  LES/RANS  models  described  above  were  used  to  simulate  turbulent  flow  over  an  Aerospatiale 
A-Airfoil  near  static  stall.  The  free- stream  Mach  number  is  0.15,  the  Reynolds  number  based  on 
a  chord  length  of  0.6m  is  2.1  x  106,  and  the  angle  of  attack  is  13.3  degrees.  An  X-Y  snapshot 
of  the  computational  mesh,  which  contains  approximately  30  x  106  cells,  is  shown  in  Figure  3.7. 
An  iso-surface  of  swirl  strength,  illustrating  the  formation  of  large  turbulent  eddies  on  the  suction 
side  of  the  airfoil,  is  illustrated  in  Figure  3.8.  Extensive  flow  measurements,  including  skin  friction, 
surface  pressure  distributions  and  velocity  and  Reynolds-stress  profiles  obtained  using  laser  Doppler 
velocimetry  measurements  are  available  for  comparison.70  Wall-resolved  large-eddy  simulations  of 
this  case  were  conducted  by  Mary  and  Sagaut.71  Other  studies  include  direct  numerical  simulations 
of  Alam  and  Sandham72  and  recent  simulations  using  a  LES/RANS  formulation  with  wall  modeling 
from  Kawai  and  Asada.73 


0  _  „  J  .  i  r  a  Figure  3.8:  Iso-surfaces  of  swirl  strength  (2000 

Figure  3.7:  X-Y  centerplane  mesh  for  A-  ...  .  .  _  _ , 

.j  s  )  illustrating  development  of  eddy  struc¬ 

tures  in  airfoil  boundary  layer. 


Figure  3.9  presents  skin- friction  and  surface  pressure  results  for  the  three  hybrid  LES/RANS  models 
tested  the  model  of  Choi,  et  al.68  which  requires  pre-calculation  of  a  model  constant,  the  model  of 
Gieseking,  et  al., 66,67  and  Gieseking’s  model  equipped  with  the  Menter-Langtry  transition  model. 
Choi’s  and  Gieseking’s  models  under-predict  the  skin  friction,  whereas  Menter-Langtry  transition 
model  under-predicts  skin-friction  on  the  front  part  of  the  airfoil  and  over-predicts  it  on  the  back 
part  of  the  airfoil.  A  large  region  of  laminar  flow  is  not  present  for  either  the  Choi  or  Gieseking 
models  but  both  show  indications  of  a  ’numerical’  transition  region.  The  use  of  the  Menter-Langtry 
transition  model  leads  to  a  large  laminar  region  terminated  by  a  laminar  separation  bubble.  Choi’s 
and  Gieseking’s  models  successfully  predict  the  pressure  plateau  near  the  trailing  edge  shown  in 
the  experimental  data,  but  Menter’s  SST  with  transition  model  fails  to  do  so.  In  the  leading  edge 
region,  the  predictions  provided  by  Gieseking’s  model  with  Menter-Langtry  are  in  better  agreement 
with  the  experimental  pressure  coefficient  data  than  those  of  Choi’s  or  Gieseking’s  models. 

Figure  3.10  shows  a  comparison  of  mean  streamwise  (tangential)  velocity  profiles  with  experimental 
measurements  at  various  axial  stations.  The  profiles  are  expressed  in  a  wall- normal  coordinate 
system.  Choi’s  model  provides  the  best  predictions  of  the  velocity  field,  with  Gieseking’s  model 
slightly  under-predicting  the  level  of  trailing-edge  separation.  The  use  of  the  Menter-Langtry  model 
promotes  flow  attachment  to  the  surface  near  the  trailing  edge. 

Streamwise  velocity  fluctuation  profiles  obtained  by  the  hybrid  LES/RANS  models  are  shown  in 
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Figure  3.9:  Mean  pressure  coefficient  (left)  and  skin  friction  coefficient  (right)  distribution  along 
the  airfoil  obtained  by  hybrid  LES/RANS  computation  compared  with  experimental  measurement. 


x/c  =  0.1  x/c  =  0.15  x/c  =  0.2  x/c  =  0.3  x/  c  =  0.5 


x/c  =  0.7  x/c  =  0.825  x/c  =  0.87  x/c  =  0.93  x/c  =  0.99 


Figure  3.10:  Mean  streamwise  velocity  profile  as  a  function  of  normalized  wall- normal  distance 
obtained  by  hybrid  LES/RANS  computation;  individual  profiles  are  separated  by  a  horizontal 
profile  of  1.4. 


Figure  3.11.  In  the  hybrid  RANS/LES  computations,  from  x/c  =  0.1  to  0.15,  the  rms  streamwise 
velocity  fluctuations  are  much  smaller  within  the  boundary  layer  when  the  transition  model  is 
included.  This  is  because  in  this  region  near  leading  edge,  the  flow  remains  laminar,  whereas  the 
flow  has  already  transitioned  and  has  become  turbulent  in  the  solutions  obtained  using  Choi’s  and 
Gieseking’s  LES/RANS  models.  From  x/c  =  0.3  to  0.5,  the  effect  of  including  the  transition  model 
is  to  enhance  the  near-surface  axial  velocity  fluctuation  intensity,  relative  to  the  other  LES/RANS 
models.  This  likely  contributes  to  the  under-prediction  of  trailing-edge  separation  mentioned  earlier. 
From  x/c  =  0.7  to  0.99,  all  models  again  under-predict  the  rms  streamwise  velocity.  Results  from 
Choi’s  model  are  in  better  agreement  with  the  experimental  data  than  those  from  Gieseking’s 
model,  whereas  Gieseking’s  model  with  Menter- Langtry  transition  leads  to  a  thinner  boundary 
layer.  This  is  not  surprising,  because  from  the  previous  mean  streamwise  velocity  analysis,  we 
know  that  the  flow  remains  attached  within  the  trailing  edge  region  when  the  transition  model  is 
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included. 


x/c  =  0.1  x/c  =  0.15  X/C  =  0.2  x/c  =  0.3  x/c  =  0.5  x/c  =  0.7  x/c  =  0.825  x/c  =  0.87  x/c  =  0.93  x/c  =  0.99 


Figure  3.11:  Profile  of  the  rms  streamwise  velocity  fluctuations  obtained  by  a  hybrid  LES/RANS 
computation;  individual  profiles  are  separated  by  a  horizontal  profile  of  0.3  (left)  and  0.2  (right). 

The  rms  wall- normal  fluctuation  velocity  profiles  obtained  by  the  hybrid  RANS/LES  models  are 
shown  in  Figure  3.12.  Trends  exhibited  in  the  LES/RANS  calculations  of  the  wall- normal  velocity 
fluctuation  generally  mirror  those  discussed  for  the  axial  velocity  fluctuation.  Fluctuation  growth 
is  suppressed,  as  expected,  near  the  leading  edge  when  the  transition  model  is  used.  Further  down¬ 
stream,  agreement  with  experiment  is  good  for  both  Choi’s  and  Gieseking’s  LES/RANS  models. 
Fluctuation  levels  in  the  near  wall  region  are  similar  among  all  models,  and  again,  the  LES/RANS 
models  provide  some  improvement  in  predictive  capability,  relative  to  the  RANS  models. 


x/c  =  0.1  x/c  =  0.15  x/c  =  0.2  x/c  =  0.3  x/c  =  0.5 


x/c  =  0.7  x/c  =  0.825  x/c  =0.87  x/c  =  0.93  x/c  =  0.99 


Figure  3.12:  Profile  of  the  rms  wall-normal  velocity  fluctuations  obtained  by  a  hybrid  LES/RANS 
computation;  individual  profiles  are  separated  by  a  horizontal  profile  of  0.3. 

Similar  trends  are  in  evidence  for  the  Reynolds  shear  stress  profiles,  as  shown  in  Figure  3.13.  Here, 
the  normalizing  factor  is  the  square  of  the  velocity,  which  tends  to  minimize  differences  among  the 
models.  Interesting,  the  use  of  the  transition  model  provides  the  best  predictions  of  the  Reynolds 
shear  stress  for  stations  up  to  x/c  =  0.5.  Further  downstream,  the  predictions  are  similar  to 
the  other  fluctuating  quantities  in  that  the  level  of  fluctuation  intensity  in  the  outer  part  of  the 
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separated  shear  layer  is  less  than  indicated  in  the  experiment.  Agreement  with  experiment  improves 
nearer  to  the  wall  for  all  models. 


x/c  =  0.1  x/c  =  0.15  x/c  =  0.2  x/c  =  0.3  x/c  =  0.5 


x/c  =  0.7  x/c  =  0.825  x/c  =0.87  x/c  =  0.93  x/c  =  0.99 


Figure  3.13:  Profile  of  the  Reynolds- averaged  shear  stress  obtained  by  hybrid  LES/RANS  compu¬ 
tation;  individual  profiles  are  separated  by  a  horizontal  profile  of  0.014. 


3.5  Interim  Conclusions 

The  results  from  the  computational  investigations  were  critical  for  the  systematic  development  of 
the  low-order  method  in  the  current  research  effort.  These  results  are  presented  in  the  following 
chapters. 
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Chapter  4 


Experimental  Facility  and  Techniques 


The  experimental  investigations  for  the  current  integrated  theoretical,  computational,  and  exper¬ 
imental  research  effort  were  performed  in  the  U.S.  Air  Force  Research  Laboratory’s  Horizontal 
Free-surface  Water  Tunnel.  This  chapter  describes  the  facility  and  the  techniques  for  force  mea¬ 
surement  and  flow  visualization. 


4.1  Facility  and  Motion  Mechanism 

The  U.S.  Air  Force  Research  Laboratory’s  Horizontal  Free-surface  Water  Tunnel  is  fitted  with  a 
three  degree  of  freedom  electric  motion  rig  enabling  independent  control  of  pitch  or  rotation,  plunge 
or  heave,  and  surge  or  streamwise- aligned  translation.  Photographs  of  the  tunnel  and  the  model 
installation  are  shown  in  Figure  4.1.  More  detail  on  the  rig  operation  is  given  in  01  et  al.8  and 
Granlund  et  ah, 74  while  the  facility  is  discussed  in  01  et  al.75 


Figure  4.1:  Test  section  and  portion  of  motion  rig  mounted  above  test  section  of  the  AFRL  Hor¬ 
izontal  Free-surface  Water  Tunnel  (left),  c  =  3”  flat  plate  with  force  balance  mounted  between 
steel  coupler  piece  and  plastic  foot  connecting  to  the  plate  (right).  Dye  injection  is  from  a  0.5  mm 
diameter  slot  at  the  plate  leading  and  trailing  edges,  |  spanwise  location. 

The  motion  rig  consists  of  a  triplet  of  H2W  linear  motors,  driven  by  AMC  DigiFlex  servo-drives 
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controlled  by  a  Galil  DMC  4040  4-channel  card,  with  user-selected  proportional/integral/derivative 
(PID)  constants  for  each  motion  channel.  Because  the  mass  of  the  linear  motors  and  plunge-rods  are 
much  larger  than  the  mass  of  test  article  and  the  expected  hydrodynamic  loads,  the  PID  constants 
are  selected  for  the  rig  to  move  itself,  and  become  largely  independent  of  the  load.  Generally 
PIDs  suitable  for  good  motion  fidelity  for  high-acceleration  motions  result  in  noise  for  smoother, 
lower-acceleration  motions,  and  as  with  all  tuning,  the  final  choice  is  a  compromise. 

The  model’s  pitch  and  plunge  are  controlled  via  two  motors  mounted  vertically  on  a  plate  above  the 
tunnel  test  section,  shown  in  the  left  portion  of  Figure  4.1.  Each  motor  actuates  a  vertical  plunge 
rod,  which  connects  via  a  bushing  to  a  coupler  piece  holding  the  force  balance.  The  upstream 
plunge  rod  is  constrained  to  move  purely  vertically,  whereas  the  downstream  plunge  rod  is  allowed 
to  pivot  in  the  test  section  vertical  plane  of  symmetry.  The  pitch  pivot  point  can  be  varied  by 
suitable  choice  of  phase  and  amplitude  difference  in  trajectory  of  front  or  rear  plunge  rod.  For 
all  cases  where  the  pitch  pivot  point  is  not  coincident  with  the  bushed  end  of  the  front  plunge 
rod,  there  will  be  a  parasitic  streamwise  displacement  of  the  plate,  which  would  be  unavoidable 
unless  the  front  plunge  rod  were  to  be  allowed  to  pivot  similarly  to  the  downstream  one.  This  is 
removed  using  the  third  degree  of  freedom,  surge,  which  also  actuates  the  fore  and  aft  translating 
motion  of  the  model.  Surge  is  achieved  using  a  larger  linear  motor  mounted  horizontally  aft  of  the 
pitch-plunge  carriage,  with  48”  peak-to-peak  stroke  and  nominal  speed  up  to  1  m/s. 

The  desired  theoretical  pitch  history  of  the  plate  is  converted  to  position  commands  for  each  linear 
motor.  Commanded  vs.  attained  displacement  histories  of  the  three  motors  are  compared  by 
interrogating  the  three  motors’  optical  encoder  tapes,  at  5000  increments/mm  for  the  two  vertical 
motors,  and  1000  increments/mm  for  the  streamwise  motor.  Peak  discrepancy  between  commanded 
and  attained  position  are  300  =t  100  increments,  which  converts  to  0.15  deg  peak  error  in  incidence 
angle,  if  the  two  vertical  motors  displacement  errors  are  in  the  worst-case  scenario  of  anti-phase. 
However,  this  is  the  incidence-angle  error  at  the  force  balance  location,  and  does  not  account  for 
vibration  or  elastic  deflection  of  the  plastic  coupler  piece  between  the  force  balance  and  the  model. 
Manual  interrogation  of  particle  image  velocimetry  raw  images  in  01  et  al.8  (not  reported  here) 
showing  reflection  of  laser  light  from  the  plate  suction-side  implies  an  upper  bound  of  incidence 
angle  uncertainty  of  0.2  deg. 


4.2  Force  Measurement 

Force  data  are  recorded  from  an  ATI  Nano-25  IP68  6-component  integral  loadcell,  oriented  with 
its  cylindrical  axis  normal  to  the  pitch-plunge-surge  plane.  It  is  capable  of  measuring  forces  in 
the  plane  of  the  airfoil  cross-section  up  to  ±125  N  and  torque  up  to  ±3  Nm.  The  published 
resolution  is  1/48  N  for  force  and  1/2640  Nm  for  torque.  The  electrical  cable  that  connects  to  the 
loadcell  is  visible  in  the  right-side  photo  in  Figure  4.1.  Load  cell  strain  gage  electrical  signals  are 
A/D  converted  in  an  ATI  NetBox  interface  and  recorded  over  Ethernet  LAN  UDP  protocol  to  a 
computer  using  a  Java  application.  The  time-base  of  the  ATI  NetBox  is  slightly  inaccurate  with 
the  clock  operating  at  a  factor  of  1.0023  too  fast  with  respect  to  physical  time.  This  is  corrected 
for  in  post-processing  of  data.  A  disadvantage  of  the  IP68  waterproofing  of  the  loadcell  is  that  it  is 
sensitive  to  immersion  depth  in  the  cylindrical  axis  direction.  Because  this  direction  is  normal  to 
the  plane  of  the  motion  of  symmetrical  models,  the  hydrostatic  force  will  not  affect  either  normal 
force,  axial  force  or  pitching  moment.  Force  and  motion  data  are  synchronized  by  polling  for  the 
trigger  signal  every  10  ms  and  starting  the  data  recording  when  initial  trigger  is  detected.  For 
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motions  that  start  with  no  position  change,  such  as  the  ramp  motions,  a  streamwise  0.001c  motion 
step  is  introduced  before  the  actual  ramp  start  in  order  to  reliably  synchronize  force  and  motion 
data. 

The  force  and  moment  signals  are  filtered  in  three  steps.  This  first  is  a  low-pass  filter  in  the 
ATI  NetBox  at  /  =  73 Hz,  to  avoid  recording  noise  not  correlated  with  motion  force  data,  but 
without  attenuating  important  fast  non-circulatory  load  spikes.  These  are  on  the  order  of  a  fifth 
as  slow  as  the  —3 dB  (half  amplitude  attenuation)  point  of  the  filter  for  the  fastest  pitch  ramp 
corner  acceleration.  The  second  step  uses  a  moving- average  of  11  points  to  smooth  the  data  while 
preserving  as  much  of  the  non-circulatory  load  spikes  as  possible.  This  smoothing  also  makes  a 
more  numerically  stable  final  step,  which  is  a  4th  order  Chebychev  II  low-pass  filter  with  —20dB 
attenuation  of  the  stopband.  The  cutoff  frequency  is  five  times  the  motion  frequency  and  the 
motion  frequency  is  calculated  by  taking  the  ramp  motion  as  a  quarter  of  a  periodic  motion.  It 
is  chosen  for  maximum  passband  flatness  and  high  rejection  of  structural  eigenfrequencies  which 
may  be  just  above  the  desired  force  frequency  information  range.  To  preclude  time-shift  of  useful 
data  in  the  passband,  the  forward-backward  filtering  technique  with  the  Matlab  filtfilt  command 
is  used. 

All  post-processing  is  done  in  Matlab.  Before  each  run,  the  loadcell  is  zero-biased  at  model  0  =  0 
deg,  which  is  adjusted  to  horizontal  with  a  bubble  level.  A  static  tare  sweep  over  0  deg<  6  <  90  deg, 
is  performed  with  500  samples  of  data  every  2  deg.  Because  the  pitch  angle  is  known  throughout 
the  motion,  and  the  position  error  is  negligible,  the  static  axial  force,  normal  force  and  pitching 
moment  due  to  static  mo  del/sting/ mount  weight  can  be  subtracted  from  the  unsteady  force  data. 


4.3  Flow  Visualization 

Flowfield  visualization  of  the  leading-edge  vortex  is  limited  to  qualitative  inferences  from  dye  in¬ 
jection.  In  planar  laser-induced  fluorescence,  a  high  concentration  of  Rhodamine  6G  in  water  is 
injected  by  a  positive-displacement  medical  infusion  pump,  connecting  to  a  set  of  0.5mm  internal- 
diameter  rigid  lines  glued  to  the  surface  of  the  plate,  spanwise  at  the  leading-  and  downstream 
at  the  trailing  edge  at  the  |  semispan  location.  The  dye  is  illuminated  by  an  Nd:YLF  527nm 
pulsed  laser  sheet  of  2mm  thickness  at  50Hz  and  images  are  recorded  with  a  PCO  DiM  x  high-ed 
camera  through  a  Nikon  PC-E  45mm  Micro  lens.  A  Tiffen  orange  #21  filter  is  used  to  remove  the 
incident  and  reflected  laser  light,  leaving  only  dye  fluorescence.  For  a  larger  spanwise  visualization, 
a  neutrally  buoyant  mixture  of  blue  food  color  and  ethanol  is  injected  and  backlight  illuminated 
by  a  Rosco  white  LED  panel.  This  method  images  the  entire  airfoil  as  well. 


4.4  Interim  Conclusions 

The  results  from  the  experimental  investigations  were  critical  for  the  systematic  development  of 
the  low-order  method  in  the  current  research  effort.  These  results  are  presented  in  Chapter  5. 
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Chapter  5 


Factors  Affecting  LEV  Initiation  in 
Unsteady  Airfoil  Flows 


In  this  chapter  we  present  a  detailed  experimental,  computational,  and  low-order  study  of  LEV 
initiation  in  airfoil  flows,  with  emphasis  on  determining  the  range  of  validity  of  the  LESP  concept. 
It  is  shown  that,  that  with  the  exception  of  slow-rate  kinematics  which  result  in  significant  trailing- 
edge  flow  separation,  the  LESP  value  for  all  other  kinematics  considered  falls  around  a  constant 
critical  value.  It  is  also  shown  that  the  LESP  concept  may  be  used  to  trigger  on-demand  or  suppress 
the  formation  of  LEVs  by  superimposing  kinematics  such  that  the  required  criteria  on  LESP  are 
met. 

5.1  Theoretical  approach 

This  section  describes  the  theoretical  methods  employed  in  this  research,  and  the  LESP  hypothesis. 
The  interested  reader  may  refer  to  references  1  and  2  for  greater  detail. 

5.1.1  Large-angle  unsteady  thin-airfoil  theory 

The  large-angle  unsteady  thin-airfoil  aims  to  eliminate  the  traditional  small-angle  assumptions  in 
thin-airfoil  theory  which  are  invalid  in  high-amplitude,  high-frequency  or  vortex-dominated  flows, 
such  as  those  considered  in  this  research.  The  theory  builds  on  the  time-stepping  approach  given 
by  Katz  V  Plotkin.76  In  Figure  5.1,  the  inertial  frame  is  given  by  OXYZ  and  the  body  frame, 
attached  to  the  moving  airfoil,  by  Bxyz.  At  time  t  =  0,  the  two  frames  coincide  and  at  time  t  >  0, 
the  body  frame  moves  toward  the  left  of  the  page  along  any  prescribed  time-varying  path  (given 
by  pitch  and  plunge  motions).  At  each  time-step,  a  discrete  trailing-edge  vortex  is  shed  from  the 
trailing  edge. 

Analogous  to  classical  thin-airfoil  theory,  the  vorticity  distribution  over  the  airfoil,  ^y(x),  is  taken 
to  be  a  Fourier  series, 


7(0)  t)  =  2[7 


.  ..IT  cos  0  v — -\  .  .  x  . 

Ao (t) — — - 1-  Y  An{t)  sm(n6) 

sin  6 

n=  1 


(1) 
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Body-fixed  Frame 
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Figure  5.1:  Illustration  of  the  time-stepping  method  in  large- angle  unsteady  thin- airfoil  theory. 


where  0  is  a  variable  of  transformation  related  to  the  chordwise  coordinate  x  as, 


x  =  -(1  —  cos#)  (2) 

Zj 

in  which  Ao(t),  Ai(t ),  ...,  An(t)  are  the  time-dependent  Fourier  coefficients,  c  is  the  airfoil  chord, 
and  U  is  the  component  of  the  airfoil’s  velocity  in  the  negative  X  direction.  The  Kutta  condition 
(zero  vorticity  at  the  trailing-edge)  is  enforced  implicitly  through  the  form  of  the  Fourier  series.  The 
Fourier  coefficients  are  determined  as  a  function  of  the  instantaneous  local  downwash  on  the  airfoil 
by  enforcing  the  boundary  condition  that  the  flow  must  remain  tangential  to  the  airfoil  surface. 


MO  =  (3) 

An(t)  =  -  f  —  cos  nOdO  (4) 

^  Jo  U 

The  induced  velocity  normal  to  the  airfoil  surface,  W(x,t),  henceforth  referred  to  as  downwash, 
is  calculated  from  components  of  motion  kinematics,  depicted  in  figure  5.2,  and  induced  velocities 
from  vortices  in  the  flowfield. 


ttt (  &V  (Tt  |  f  •  .  d&tev  \  tt  •  •/  \  |  f  d&tev 

W (x,  t)  =  — —  =  —  (U  cos  a  +  n  sm a  H - - — )  —  U  since  —  a(x  —  ac)  +  ri  cos  a - - —  (5) 

oz  ox  ox  oz 

where  <j)B  and  (j)tev  are  the  velocity  potentials  associated  with  bound  and  trailing-edge  vorticity, 
rj(x)  is  the  camber  distribution  on  the  airfoil,  d^v  and  are  velocities  induced  tangential  and 
normal  to  the  chord  by  trailing  edge  discrete  vortices.  The  motion  parameters  include  the  plunge 
velocity  in  the  Z  direction,  /i,  and  the  pitch  angle  of  the  chord  with  respect  to  the  X  direction, 
a.  Trailing-edge  vortices  are  shed  at  every  time-step  as  mentioned  earlier,  and  their  strengths  are 
calculated  iteratively  such  that  Kelvin’s  circulation  condition  is  enforced. 
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Figure  5.2:  Airfoil  velocities  (positive  as  shown)  and  pivot  location. 


Ntev 

^  r tevm  =  o 

m=  1 


(6) 


where  is  the  bound  circulation  calculated  by  integrating  the  chordwise  distribution  of  bound 
vorticity  over  the  airfoil  chord: 


Tb  =  Uctt 


Ao(t)  + 


Mjt) 

2 


(7) 


5.1.2  Leading  Edge  Suction  Parameter  (LESP) 


Figure  5.3:  Depiction  of  flow  around  a  thin  airfoil’s  leading  edge. 

In  thin-airfoil  theory,  the  airfoil  thickness  and  hence  the  leading-edge  radius  is  zero.  This  requires 
the  flow  to  turn  180°  around  the  leading  edge  (figure  5.3),  giving  rise  to  a  theoretically  infinite  flow 
velocity  at  the  leading  edge,  Vle,  of  a  thin  airfoil.  From  Garrick77  and  von  Karman  &  Burgers,78 
we  have  that  the  form  of  this  theoretically  infinite  velocity  is  given  by, 

=  lim 

x — >LE  X 

where  S  is  a  measure  of  the  suction  at  the  leading  edge  and  is  given  by, 


(8) 


23 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


(9) 


S  =  lim  -7 (x,t)\/x 

x^LE  2 

Since  7 (x,t)  is  infinite  in  order  of  l/^/x  at  the  leading  edge,  the  value  of  S  is  finite.  Evaluating 
using  the  current  formulation, 

S  =  y/cUAo(t)  (10) 

Because  the  LESP  is  a  nondimensional  measure  of  the  suction  at  the  leading  edge,  5,  for  given 
values  of  c  and  U,  we  may  simply  equate  it  to  the  Ao  value  as, 

LESP(t)  =  A0(t)  (11) 

As  discussed  in  Katz,20  real  airfoils  have  rounded  leading  edges  which  can  support  some  suction 
even  when  the  stagnation  point  is  away  from  the  leading  edge.  The  amount  of  suction  that  can  be 
supported  is  dependent  on  the  airfoil  shape  and  Reynolds  number  of  operation.  Since  the  LESP 
(the  Ao  value)  is  a  measure  of  the  suet  ion/ velocity  at  the  leading  edge,  it  is  a  logical  choice  to 
develop  a  correlation  for  initiation  of  LEV  formation  based  on  the  LESP. 

Using  an  inviscid  parameter  to  predict  trends  in  viscous  behavior  is  not  a  new  idea.  For  example, 
Rival  et  al.13  have  shown  that  LEV  detachment/”  pinch-off”  occurrs  when  the  stagnation  point  of 
the  LEV  breaches  the  airfoil’s  trailing  edge  -  a  phenomenon  that  may  be  predicted  from  inviscid 
theory.  The  Ao  term  has  also  been  previously  used  to  develop  useful  correlations  in  steady  aerody¬ 
namic  theory.  It  is  well  known,  for  instance,  that  the  ideal  lift  coefficient  of  a  laminar-flow  airfoil 
in  steady  flow,  which  usually  falls  close  to  the  middle  of  the  drag  bucket,  corresponds  to  the  lift 
coefficient  at  which  the  Ao  coefficient  is  zero;79,80  this  idea  can  be  used  to  estimate  the  C/ -shift  in 
the  drag  bucket  due  to  a  trailing-edge  cruise  flap.81  The  Ao  term  in  thin-airfoil  theory  is  the  only 
term  that  results  in  a  singularity  in  the  vorticity  distribution  at  the  leading  edge,  and  hence  is  a 
good  measure  of  the  flow  at  the  leading  edge. 

The  research  of  Morris  V  Rusak82  on  inception  of  leading-edge  stall  on  stationary,  two-dimensional, 
smooth,  thin  airfoils  provides  another  explanation  of  why  the  critical  value  of  the  Ao  term  should 
correspond  to  initiation  of  LEV  formation.  In  their  work,  the  authors  have  used  matched  asymptotic 
theory  with  the  flow  around  most  of  the  airfoil  chord  described  in  terms  of  an  outer  region  which 
is  solved  using  thin  airfoil  theory.  The  flow  in  the  vicinity  of  the  leading  edge  forms  the  inner 
region,  which  is  treated  as  a  model  problem  of  a  uniform,  incompressible  and  viscous  flow  past 
a  semi-infinite  parabola  and  solved  through  numerical  simulations  of  the  unsteady  Navier-Stokes 
equations.  The  flows  in  the  inner  and  outer  regions  are  made  to  asymptotically  match  each  other. 
The  far-field  circulation  for  the  inner  flow  is  governed  by  a  parameter  that  is  related  to  the  airfoil’s 
angle  of  attack.  This  approach  allows  the  determination  of  the  critical  angle  of  attack  for  leading- 
edge  stall  onset  as  the  condition  at  which  a  global  separation  zone  is  predicted  in  the  solution  for 
the  inner  flow.  For  a  given  airfoil  geometry,  the  Ao  value  is  linearly  related  to  the  angle  of  attack 
in  stationary  flow.  This  shows  that,  at  a  given  Reynolds  number,  leading-edge  stall  in  stationary 
airfoil  flow  is  related  to  a  critical  value  of  the  Ao  coefficient.  In  the  current  unsteady  thin  airfoil 
theory,  the  Ao  value  accounts  not  only  for  the  instantaneous  angle  of  attack,  but  also  for  the  motion 
kinematics  and  the  effect  of  vorticity  in  the  flow  through  the  zero-normal  flow  boundary  conditions 
in  equation  5.  It  follows,  therefore,  that  the  critical  value  of  Ao  would  correspond  to  initiation  of 
LEV  formation  in  unsteady  flow. 
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5.2  Study  of  LEV  initiation  in  low  Reynolds-number  flows 


Trends  in  LEV  initiation  for  various  motion  kinematics  and  the  validity  of  the  LESP  hypothesis 
are  analyzed  in  this  section  using  data  sets  from  CFD  and  experiments.  An  SD7003  airfoil  at  a 
Reynolds  number  of  20,  000  is  used  for  all  cases,  and  a  parameter  space  encompassing  a  range  of 
values  for  various  kinematic  terms  is  chosen  as  described  in  section  5.2.1. 


5.2.1  Definition  of  motion  kinematics 


Combinations  of  pitch  and  plunge  maneuvers  are  considered  for  the  parametric  study.  Both  these 
motions  are  generated  with  a  modified  version  of  the  Eldredge  function  which  produces  a  ramp 
motion  with  smoothed  corners.  The  pitch  histories  are  given  by: 


K* 

—  &start  V 

CLs 


cosh  (as(t*  —  £*)) 
cosh  (as(t*  —  £2)) 


a 


amp 


(12) 


,  where  as  is  the  smoothing  parameter  defined  as: 


CLS 


n2Kr 


1  cr) 


(13) 


and, 


£2  =  *1  + 


a 


amp 


2  Kr 


(14) 


In  these  equations,  t\  denotes  the  time  at  start  of  ramp,  and  denotes  the  time  at  end  of  ramp. 
In  all  simulations  in  this  paper,  £*  is  taken  as  5.0  to  generate  a  steady  starting  solution  and  hence 
minimize  the  effect  of  starting  vortices  on  the  solution.  The  parameter  a  is  a  nondimensional 
measure  of  smoothing,  and  is  equal  to  0.8  in  all  kinematics  considered  here.  Ka  is  the  reduced 
frequency  of  pitch.  The  term  astart  is  used  to  generate  kinematics  where  the  ramp  starts  from  a 
non-zero  value. 

The  plunge  kinematics  are  constructed  with  the  same  equations,  by  replacing  a  with  h/c,  Ka  with 
Kh  and  aamp  with  (h/c)amp.  hstart  is  not  used  (always  0).  As  the  plunge  motion  is  in  combination 
with  pitch  in  this  study,  the  reduced  frequency  for  plunge  is  chosen  such  that  the  pitch  and  plunge 
ramps  occupy  the  same  nondimensional  time  (Kh  =  Ka  *  ( h/c)arnp/aarnp ).  Hence,  equations  13 
and  14  are  not  altered.  The  variation  in  plunge  is  given  by, 


h  Ka{h/  c)amp  cosh(as(£  £1)) 
c  asaamp  [cosh(as(£*  -  £2)) 

A  typical  pitch-only  ramp  motion,  starting  at  astart  =  0,  with  pitch  amplitude  aamp  =  30  deg, 
reduced  frequency  Ka  =  0.2  and  pivoted  at  quarter-chord  is  considered  as  a  baseline  for  the 
parametric  study.  The  detailed  definition  of  this  case  is  given  in  table  5.1. 

The  parameter  space  for  case  studies  1-3  is  constructed  by  considering  variations  in  pitch-axis 
location,  pitch  reduced  frequency  and  pitch  start  angles  with  respect  to  the  baseline  case.  In  case 
study  4,  a  pitch-plunge  combination  is  compared  to  the  baseline  case.  In  keeping  with  the  bounds 
of  the  experimental  apparatus  with  respect  to  plunge  maneuvers,  the  Ka  for  the  pitch-plunge 


+ 


Wc) 


amp 
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Table  5.1:  Base  parameter  set  used  to  study  LEV  initiation  in  low- Reynolds  number  flow 


Parameter  Symbol  Value 


Reynolds  number 
Start  time  of  ramp 
Smoothing  parameter 
Initial  pitch  angle 
Pitch  amplitude 
Pitch  rate 
Pivot  location 
Plunge  amplitude 
Plunge  rate 


Re 

20,000 

li 

5.0 

cr 

0.8 

start  (deg) 

0 

&  amp  (deg) 

30 

II 

§18' 

0.20 

Xp/x 

0.25 

(h/ c)amp 

0.0 

II 

0.0 

combination  was  chosen  to  be  0.1  rather  than  the  baseline  value  of  0.2.  As  mentioned  previously, 
the  Kh  is  calculated  such  that  the  pitch  and  plunge  ramps  occupy  the  same  nondimensional  time 
(Kh  =  Ka  *  (h/ c) arnp/ a arnp).  The  kinematics  for  the  4  studies  are  listed  in  table  5.2. 

Table  5.2:  Parameter  variations  used  to  study  LEV  initiation  in  low- Reynolds  number  flow  (baseline 
values  are  in  bold) 


Case  study 

Variable  Parameter 

Values 

1 

Pivot  location  (xp/c) 

0.0,  0.25,  0.75 

2 

Pitch  rate  (Ka) 

0.01,  0.03,  0.05,  0.1,  0.2,  0.4 

3 

Initial  pitch  angle  ( a  start ) 

10,  5,  0,  -5,  -10, -15 

4 

Pitch-plunge  combination 

OLamp  30,  Ka  0.1,  (h/c)amp  — 0.1,  Kh  0.0191 

For  all  kinematics,  the  pitch  angles  at  which  LEV  formation  is  initiated  are  determined  using 
experiments  and  CFD,  and  the  LESP  values  at  these  pitch  angles  are  determined  from  theory. 
Skin-friction  coefficient  data  from  CFD  is  used  to  quantitatively  identify  the  initiation  of  LEV 
formation.  Experimental  data  is  used  to  mutually  validate  the  CFD  and  qualitatively  study  the 
initiation  of  LEV  formation  for  all  kinematics.  The  identification  of  LEV  initiation  using  skin- 
friction  data  is  illustrated  with  the  baseline  case  in  section  3.2.  Four  parametric  studies  are  then 
considered  in  section  5.3  to  study  LEV  initiation  behavior  and  to  establish  the  envelope  of  validity 
of  the  LESP  hypothesis.  Section  5.3.5  presents  the  results  from  all  4  case  studies  in  combination. 
Finally,  section  5.3.6  discusses  two  design  cases  in  which  plunge  motions  are  added  to  the  baseline 
pitch  motion  to  alter  the  occurrence  of  LEV  formation. 


5.3  Results  of  parametric  studies  of  LEV  initiation 

5.3.1  Case  study  1:  kinematics  with  varying  pivot  locations 

In  this  study,  the  effect  of  pitch-axis  location  on  LEV  initiation  in  unsteady  maneuvers  is  examined. 
In  addition  to  the  baseline  case  (quarter-chord  pivot),  two  cases  with  modified  pivot  locations 
(leading-edge  and  three-quarter-chord)  are  considered.  The  pitch-angle  during  the  ramp-maneuver 
at  which  LEV  formation  is  initiated  (as  determined  using  the  procedure  described  in  section  3.2) 
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Figure  5.4:  Case  study  1:  pitch-angle  varia¬ 
tion  with  time,  and  pitch-angles  at  which  LEV 
formation  is  initiated  (kinematics  with  different 
pivot  locations). 


Figure  5.5:  Case  study  1:  LESP  variation  with 
time,  and  LESP  values  at  which  LEV  forma¬ 
tion  is  initiated  (kinematics  with  different  pivot 
locations). 


is  plotted  for  all  three  cases  in  figure  5.4.  As  the  pivot  location  is  varied  along  the  chord  from 
leading  to  trailing  edge,  the  pitch  angle  at  LEV  initiation  is  seen  to  increase.  Clearly,  there  is  no 
obvious  relation  between  initiation  of  LEV  formation  and  the  values  of  the  pitch  angle  at  that  time 
instant.  In  figure  5.5,  the  time  variation  of  LESP  for  the  3  kinematics  from  unsteady  thin-airfoil 
theory  are  co-plotted  with  the  instants  of  LEV  formation  marked.  It  is  seen  that  the  initiation 
of  upper-surface  LEV  formation  occurs  at  a  near-constant  LESP  value.  The  results  for  this  case 
study  hence  confirm  the  LESP  hypothesis  that  rounded  edges  can  support  a  certain  maximum 
suction  force,  beyond  which  LEV  formation  occurs. 

More  insight  into  the  LEV  formation  and  shedding  process  is  gained  by  analyzing  the  experimen¬ 
tal  and  numerical  results  together.  Flow  visualization  plots  from  experiment,  vorticity  plots  from 
CFD,  pressure  coefficient  (upper  and  lower  surfaces)  and  skin- friction  coefficient  (upper  surface) 
distributions  from  CFD  at  the  instants  of  LEV  initiation  for  the  three  cases  are  shown  in  figure  5.6. 
At  the  time  of  LEV  initiation,  we  see  that  the  Cp  and  Cf  distributions  for  the  3  cases  are  qualita¬ 
tively  similar.  The  vorticity  and  flow  visualization  plots  show  that  the  flow  is  attached  over  most 
of  the  airfoil  for  all  3  cases,,  as  is  expected  for  the  high  reduced  frequency  K  =  0.2.  This  explains 
the  LESP  hypothesis  holding  valid  in  this  case  study,  as  the  theory  is  derived  on  the  basis  of  an 
attached  flow  assumption. 

The  same  plots  from  CFD  and  experiment,  at  some  time  after  the  initiation  of  LEV  formation  are 
shown  in  figure  5.7  .  The  time  instant  is  chosen  as  that  when  the  inviscid  LESP  value  is  greater 
than  LESPcrit  by  a  value  of  0.1.  The  LESP  is  a  measure  of  the  velocity  at  the  leading  edge  and 
also  the  velocity  at  the  start  of  the  shear  later  emanating  from  the  leading  edge.  Now,  the  flux 
of  vorticity  into  the  free  shear  layer  (LEV)  may  be  related  to  the  shear  layer  edge  velocity  by  the 
experimentally  observed  relation  of  Fage  and  Johansen,83 

dL Sh  1  t  y2 

~dT  “  2  sh 
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Figure  5.6:  Case  study  1:  left  to  right  -  flow  visualization  from  experiment,  vorticity  plots  from 
CFD  ,  Cp  (upper  and  lower  surfaces)  and  Cf  (upper  surface)  distributions  from  CFD  at  the  instants 
of  LEV  initiation.  Top  to  bottom  -  pivot  location  at  leading  edge,  quarter  chord  and  three-quarter 
chord. 


This  relation  has  also  been  used  to  model  vortex  shedding  from  edges  in  discrete-vortex  methods 
such  as  those  by  Sarpkaya17  and  Katz.20  The  instants  when  the  LESP  is  greater  than  LESPcrit 
by  a  certain  constant  value  thus  correspond  to  times  when  approximately  equal  vorticity  has  been 
shed  into  the  LEVs  in  all  cases.  Strictly,  the  LESP  curves  in  figure  5.5  are  only  valid  until  the 
instant  of  LEV  initiation,  as  LEV  shedding  which  occurs  after  this  instant  is  not  modeled.  Still, 
A  LESP  =  0.1  is  chosen  as  a  measure  when  all  cases  would  be  at  approximately  the  same  stage  of 
vortex  development,  irrespective  of  motion  parameters.  Figure  5.7  shows  that  the  LEV  structures 
and  Cp  distributions  are  similar  for  the  three  cases  with  different  pitch-axis  locations.  Further,  all 
3  cases  show  the  presence  of  a  concentrated  LEV  with  no  noticeable  flow  separation  over  the  rest 
of  the  airfoil. 

5.3.2  Case  study  2:  kinematics  with  varying  pitch  rates 

In  this  case  study,  the  effect  of  pitch-rate  on  LEV  formation  in  unsteady  maneuvers  is  investigated. 
In  addition  to  the  baseline  case  (K  =  0.2),  five  cases  with  modified  pitch  rates  (0.01,  0.03,  0.05, 
0.1  and  0.4)  are  considered. 

The  pitch-angle  during  the  ramp  motion  at  which  LEV  formation  is  initiated  (as  determined  using 
the  procedure  described  in  section  3.2)  is  plotted  for  all  six  cases  in  figure  5.8.  The  pitch  angle  at 
LEV  initiation  is  seen  to  increase  with  increasing  pitch  rate.  High-pitch  rates  hence  serve  to  keep 
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Figure  5.7:  Case  study  1:  left  to  right  -  flow  visualization  from  experiment,  vorticity  plots  from  CFD 
,  Cp  (upper  and  lower  surfaces)  and  Cf  (upper  surface)  distributions  from  CFD  at  A LESP  =  0.1 
after  the  instants  of  LEV  initiation.  Top  to  bottom  -  pivot  location  at  leading  edge,  quarter  chord 
and  three-quarter  chord. 


the  flow  attached  to  the  airfoil,  which  is  well  known  from  dynamic  stall  research  (for  e.g.  ref84).  In 
figure  5.9,  the  time  variation  of  LESP  for  the  six  kinematics  from  unsteady  thin- airfoil  theory  are 
co-plotted  with  the  instants  of  LEV  formation  marked.  We  see  that  initiation  of  LEV  formation 
for  cases  with  relatively  high  pitch  rates  (K  >=  0.1)  occurs  at  a  near  constant  LESP  value.  The 
LESP  values  at  LEV  initiation  for  the  cases  with  lower  pitch  rates  (K  <  0.1)  are  lower,  with  the 
values  decreasing  for  lower  pitch  rates.  This  result  is  analyzed  in  detail  by  studying  flow  features 
from  experiments  and  CFD  below. 

Flow  visualization  plots  from  experiment,  vorticity  plots  from  CFD,  pressure  coefficient  (upper  and 
lower  surfaces)  and  skin-friction  coefficient  (upper  surface)  distributions  from  CFD  at  the  instants 
of  LEV  initiation  for  the  six  cases  are  shown  in  figure  5.10.  The  first  three  cases  with  pitch  rates 
of  0.01,  0.03  and  0.05  respectively,  are  seen  to  exhibit  significant  boundary- layer  thickening  and 
flow  separation  on  the  airfoil  surface  at  the  time  of  LEV  initiation.  For  these  three  cases,  the  flow 
visualization  from  experiment  shows  that  flow  is  separated  over  more  than  50%  of  the  airfoil  at 
the  time  of  LEV  initiation.  Hence  we  do  not  expect  the  LESP  hypothesis  to  hold  true  for  these 
cases,  as  the  underlying  unsteady  thin-airfoil  theory  assumes  attached  flow  over  the  whole  airfoil 
and  is  hence  not  valid.  If  trailing-edge  separation  were  modeled  in  the  calculation  of  LESP ,  the 
hypothesis  may  hold  true  even  for  the  slow  pitch-rate  cases.  The  Cp  distributions  for  the  six  cases 
show  a  clear  trend,  with  higher  pitch  rates  resulting  in  a  greater  values  of  pressure  coefficient  on  the 
airfoil.  The  Cf  distributions  and  the  Cy-spike  which  is  used  to  identify  LEV  initiation  also  show 
clear  trends  with  the  spike  moving  aft  on  the  airfoil  and  getting  broader /diffused  with  decreasing 
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Figure  5.8:  Case  study  2:  pit  eh- angle  varia¬ 
tion  with  time,  and  pitch-angles  at  which  LEV 
formation  is  initiated  (kinematics  with  different 
pitch  rates). 


Figure  5.9:  Case  study  2:  LESP  variation  with 
time,  and  LESP  values  at  which  LEV  forma¬ 
tion  is  initiated  (kinematics  with  different  pitch- 
rates)  . 


pitch  rate. 

The  same  plots  from  CFD  and  experiment,  at  some  time  after  the  initiation  of  LEV  formation 
are  shown  in  figure  5.7.  The  time  instant,  as  earlier,  is  chosen  as  that  when  the  inviscid  LESP 
value  is  greater  than  LESPcrn  by  a  value  of  0.1.  The  LEV  structures  for  the  six  cases  exhibit  a 
trend  of  being  more  concentrated  with  increasing  pitch  rate.  For  the  first  case  with  K  =  0.01,  the 
pitch-rate  is  so  slow  that  the  airfoil  is  “nearly  steady”  and  displays  a  bluff-body  type  flow.  The 
second  and  third  cases  show  a  distinct  LEV,  along  with  the  presence  of  flow  separation  over  the 
rest  of  the  airfoil.  The  final  three  cases  have  a  concentrated  LEV  with  no  significant  separation 
over  the  rest  of  the  airfoil  surface.  From  figure  5.9,  we  see  that  these  cases  display  LEV  initiation 
at  a  near  constant  LESP  value,  whereas  the  LESP  at  LEV  initiation  for  the  slower  cases  occurs  at 
progressively  lower  values. 


5.3.3  Case  study  3:  kinematics  with  varying  initial  pitch  angles 

In  case  study  3,  the  influence  of  boundary  layer  development  and  flow  separation  on  LEV  initiation 
is  investigated  by  starting  the  pitch  ramps  from  non-zero  values.  In  addition  to  the  baseline  case 
(which  starts  from  a  pitch  angle  of  0),  five  cases  starting  with  starting  pitch  angles  of  10,  5,  —5, 
—  10  and  —15  deg  are  used,  from  The  pitch  angle  histories  for  these  cases  and  the  angles  at  which 
LEV  formation  is  initiated  are  plotted  for  all  six  cases  in  figure  5.12. 

In  all  cases,  a  steady  boundary  layer  is  established  at  the  starting  pitch  angle  before  the  ramp  is 
initiated.  The  results  show  some  variation  in  the  pitch  angles  at  which  LEV  formation  is  initiated 
for  these  cases.  There  is  no  clear  trend,  suggesting  that  the  influence  of  an  established  boundary 
layer  on  LEV  formation  is  nonlinear.  Figure  5.13  displays  the  time  variation  of  LESP  as  determined 
from  theory  for  the  6  cases.  For  cases  starting  with  non-zero  pitch  angles,  the  LESP  values  at  LEV 
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Figure  5.10:  Case  study  2:  left  to  right  -  flow  visualization  from  experiment,  vorticity  plots  from 
CFD  ,  Cp  (upper  and  lower  surfaces)  and  Cf  (upper  surface)  distributions  from  CFD  at  the  instants 
of  LEV  initiation.  Top  to  bottom  -  pitch  rate  of  0.01,  0.03,  0.05,  0.1,  0.2  and  0.4). 


initiation  are  lower  than  that  for  the  baseline  case  (similar  to  trend  seen  in  section  5.3.2,  though 
they  are  still  close  to  the  baseline  value. 

Flow  visualization  plots  from  experiment,  vorticity  plots  from  CFD,  Cp  (upper  and  lower  surfaces) 
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Figure  5.11:  Case  study  2:  left  to  right  -  flow  visualization  from  experiment,  vorticity  plots  from 
CFD  ,  Cp  (upper  and  lower  surfaces)  and  Cf  (upper  surface)  distributions  from  CFD  at  A LESP  = 
0.1  after  the  instants  of  LEV  initiation.  Top  to  bottom  -  pitch  rate  of  0.01,  0.03,  0.05,  0.1,  0.2  and 
0.4). 


and  Cf  (upper  surface)  distributions  from  CFD  at  the  instants  of  LEV  initiation  for  cases  study 
3  are  shown  in  figure  5.14.  The  first  2  cases  (starting  at  10  and  5  degrees)  exhibit  flow  separation 
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Figure  5.12:  Case  study  3:  pitch- angle  varia¬ 
tion  with  time,  and  pitch-angles  at  which  LEV 
formation  is  initiated  (kinematics  with  different 
start  angles). 


Figure  5.13:  Case  study  3:  LESP  variation  with 
time,  and  LESP  values  at  which  LEV  formation 
is  initiated  (kinematics  with  different  start  an- 
gles). 


on  the  airfoil  upper  surface  at  the  time  of  LEV  initiation,  while  the  final  3  cases  (starting  at  —5, 
—  10  and  —15  deg)  exhibit  separated  flow  on  the  lower  surface.  The  LESP  values  at  the  instants  of 
LEV  formation  for  these  cases  were  however  quite  close  as  observed  in  figure  5.13. 

The  same  plots  from  CFD  and  experiment,  at  an  instant  after  the  initiation  of  LEV  formation 
when  the  inviscid  LESP  value  is  greater  than  LESPcra  by  a  value  of  0.1,  are  shown  in  figure  5.7. 
All  the  cases  exhibit  a  concentrated  vortex,  with  the  first  two  also  showing  significant  trailing-edge 
separation.  The  LESP  hypothesis  is  however  seen  to  hold  for  all  cases  (figure  5.13)  in  contrast  to 
the  slow  ramp  cases  in  section  5.3.2  which  also  had  significant  trailing-edge  separation.  Hence  the 
presence  of  separated  shear  layers/thick  boundary  layers  as  an  initial  condition  does  not  appear  to 
affect  LEV  initiation  as  dynamics  flow  separation  on  the  airfoil  surface  progressing  from  the  trailing 
edge  towards  the  leading  edge  owing  to  low  pitch  rates. 


5.3.4  Case  study  4:  kinematics  with  pitch-plunge  combination 

In  this  study,  we  aim  to  establish  that  the  LESP  hypothesis  (LEV  initiation  occurring  at  the 
same  critical  value  of  LESP)  applies  not  only  to  various  pitching  maneuvers,  but  to  any  arbitrary 
unsteady  maneuver.  In  addition  to  the  baseline  case,  a  pitch-plunge  combination  is  considered.  The 
latter  has  a  pitch  amplitude  of  30  deg,  plunge  amplitude/ chord  of  —0.1,  and  reduced  frequency  (in 
both  pitch  and  plunge)  of  0.1  as  given  in  table  5.2.  The  pitch  angle  histories  for  the  two  cases  and 
the  angles  at  which  LEV  formation  is  initiated  are  plotted  in  figure  5.16. 

Figure  5.17  displays  the  time  variation  of  LESP  as  determined  from  theory  for  the  two  cases, 
and  the  LESP  values  at  the  instant  of  LEV  initiation  as  determined  from  CFD  (section  3.2.  LEV 
initiation  in  both  cases  is  seen  to  occur  at  the  same  LESP  value,  thereby  validating  the  LESP 
hypothesis  for  arbitrary  kinematics. 
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Figure  5.14:  Case  study  3:  left  to  right  -  flow  visualization  from  experiment,  vorticity  plots  from 
CFD  ,  Cp  (upper  and  lower  surfaces)  and  Cf  (upper  surface)  distributions  from  CFD  at  the  instants 
of  LEV  initiation.  Top  to  bottom  -  starting  pitch  angles  of  10,  5,  0,  —5,  —10  and  —15  deg). 
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Figure  5.15:  Case  study  3:  left  to  right  -  flow  visualization  from  experiment,  vorticity  plots  from 
CFD  ,  Cp  (upper  and  lower  surfaces)  and  Cf  (upper  surface)  distributions  from  CFD  at  A LESP  = 
0.1  after  the  instants  of  LEV  initiation.  Top  to  bottom  -  starting  pitch  angles  of  10,  5,  0,  —5,  —10 
and  —15  deg). 
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Figure  5.16:  Case  study  4:  pitch  angle/plunge 
amplitude  variation  with  time,  and  pitch  angles 
at  which  LEV  formation  is  initiated  (baseline 
case  and  pitch-plunge  combination). 


Figure  5.17:  Case  study  4:  LESP  variation  with 
time,  and  LESP  values  at  which  LEV  formation 
is  initiated  (baseline  case  and  pitch-plunge  com¬ 
bination). 


As  for  the  other  case  studies,  flow  visualization  from  experiment,  vorticity  plots  from  CFD,  Cp, 
and  Cf  distributions  from  CFD  at  the  instants  of  LEV  initiation  for  the  two  cases  are  shown  in 
figure  5.18.  Despite  LEV  initiation  occurring  at  different  angles  of  attack,  the  Cf  distributions 
for  the  two  cases  are  similar  at  the  instant  of  LEV  initiation,  which  also  corresponds  to  the  same 
critical  value  of  LESP  as  seen  in  figure  5.17. 

The  same  plots  from  CFD  and  experiment  for  the  two  cases,  at  a  time  after  initiation  of  LEV 
formation  when  the  inviscid  LESP  value  is  greater  than  LESPcrit  by  a  value  of  0.1,  are  shown  in 
figure  5.19.  The  vortex  development  for  the  two  cases  is  seen  to  be  quite  different  owing  to  the 
different  types  of  motion  and  the  different  reduced  frequencies.  While  the  baseline  case  exhibits  a 
small  concentrated  vortex,  the  pitch-plunge  combination  evinces  a  more  diffused  LEV  accompanied 
by  tr ailing-edge  flow  separation. 

5.3.5  Compilation  of  all  test  cases 

The  results  from  all  parametric  studies  performed  in  this  research  are  compiled  here.  The  applica¬ 
bility  of  the  LESP  hypothesis  for  predicting  LEV  initiation  is  demonstrated  using  figure  5.20.  On 
the  y-axis,  are  the  predictions  from  CFD  for  angles  of  attack  at  LEV  initiation  in  all  the  kinematics 
considered.  On  the  x-axis,  are  angles  of  attack  predicted  from  inviscid  theory,  using  the  critical 
LESP  value.  The  LESP  value  of  the  baseline  case  at  the  instant  of  the  LEV  initiation  is  taken  to 
be  the  critical  LESP.  The  instant  of  LEV  initiation  for  any  other  case  is  predicted  from  theory  as 
the  instant  when  the  instantaneous  LESP  value  just  crosses  the  critical  LESP  value.  As  seen  from 
figure  5.20,  with  the  exceptions  of  the  three  slow  pitch-rate  cases  from  case  study  2  (K  =  0.01, 
0.03  and  0.05),  the  predictions  from  CFD  and  theory  are  within  an  error  margin  of  ±1  deg  for 
most  cases,  and  within  ±2.5  deg  for  all  cases.  This  clearly  demonstrates  the  viability  of  the  LESP 
concept  as  a  tool  for  predicting  LEV  initiation  using  inviscid  theory.  The  critical  LESP  needs  to 
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Figure  5.18:  Flow  visualization  from  experiment,  vorticity  plots  from  CFD  ,  Cp  (upper  and  lower 
surfaces)  and  Cf  (upper  surface)  distributions  from  CFD  at  the  instant  of  LEV  initiation  for  cases 
1  and  15. 


Figure  5.19:  Flow  visualization  from  experiment,  vorticity  plots  from  CFD,  Cp  (upper  and  lower 
surfaces)  and  Cf  (upper  surface)  distributions  from  CFD  at  some  A LESP  after  the  initiation  of 
LEV  formation  for  cases  1  and  15. 


37 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Figure  5.20:  Comparison  of  angle  of  attack  for  LEV  initiation  from  low-order  and  high-order 
methods. 


be  calibrating  using  just  one  motion  for  a  given  airfoil  shape  and  Reynolds  number  combination. 
It  can  then  be  employed  to  predict  LEV  initiation  for  any  other  (fast-rate)  motion  with  the  same 
airfoil  and  Reynolds  number  combination. 

Figure  5.21  shows  the  LESP  values  when  LEV  initiation  occurs  for  all  the  case  studies  considered, 
and  further  reinforces  the  LESP  hypothesis.  With  the  exception  of  the  3  slow  pitch-rate  cases  in 
case  study  2  (K  =  0.01,  0.03  and  0.05),  the  LESP  values  for  all  other  cases  at  the  instant  of  LEV 
initiation  are  seen  to  occur  around  the  same  critical  value.  The  figure  shows  that  with  the  exception 
of  the  three  outliers,  the  LESP  values  at  LEV  initiation  for  all  cases  lies  within  14%  of  the  critical 
value.  In  the  next  section,  the  use  of  LESP  as  a  low-order  tool  for  manipulating  LEV  formation 
is  demonstrated. 

5.3.6  Trigger /Suppress  LEV  formation  using  LESPcrit 

In  the  preceding  sections,  we  showed  that  the  start  of  separation  at  the  leading  edge  is  related  to  the 
LESP  exceeding  a  certain  critical  value.  The  critical  LESP  value  is  a  function  of  the  airfoil  shape 
and  Reynolds  number  of  operation.  Once  pre-determined  using  experimental  or  computational 
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Figure  5.21:  LESP  values  at  the  instant  of  LEV  initiation  (as  determined  from  CFD),  compiled  for 
all  cases  considered  . 


methods,  the  critical  LESP  value  corresponds  to  onset  of  separation  at  the  leading  edge,  irrespective 
of  motion  kinematics.  Hence  LEV  occurrence  may  be  controlled  by  suitably  altering  the  motion 
kinematics  such  that  the  LESP  critical  value  is  attained  at  the  required  time  of  start  of  LEV 
formation. 

Consider  the  baseline  case  in  section  3.2,  analyzed  with  CFD  -  a  pitch  motion  with  amplitude  of 
30  deg,  reduced  frequency  of  0.2,  and  pivot  about  the  quarter  chord.  For  these  motion  kinematics, 
LEV  formation  was  seen  to  be  initiated  on  the  upper  surface  at  t*  =  0.95,  a  =  22  deg.  In  this 
section,  we  modify  the  occurrence  of  LESP  variation  by  superimposing  a  plunge  motion  on  the 
baseline  case,  such  that  LEV  initiation  is  either  advanced  or  delayed  as  desired. 

As  it  is  well  known  that  (negative)  rate  of  plunge  is  equivalent  to  variation  in  pitch  angle,7  the 
plunge  rate  is  taken  to  be  in  the  form  of  an  Eldredge-ramp  function  (same  as  pitch  angle). 

h  _  Ka(h/c)amp  "cosh (a3(t*  -  £*)) 
c  asaamp  [cosh(as(t*  - 1^)) 

The  plunge  motion  to  be  superimposed  on  pitch  is  obtained  by  integrating  the  Eldredge-form  plunge 
rate, 


+ 


[h/cJ 


amp 


(17) 


h 

c 


(18) 


A  Newton  iteration  is  used  to  determine  the  value  of  (h/c)  such  that  criteria  on  LESP  are  satisfied, 
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Figure  5.22:  Pitch  amplitude  variation  with 
time  for  the  baseline  case,  and  plunge  ampli¬ 
tude  variations  which  are  used  in  combination 
with  baseline  pitch  to  generate  the  two  design 
cases. 


Figure  5.23:  LESP  variation  with  time  for  the 
baseline  case  and  the  two  design  cases.  As  re¬ 
quired  by  the  design  criteria,  the  intersection  of 
instantaneous  LESP  with  critical  LESP  is  at 
t*  —  t\  =  0.5  and  1.5  respectively. 


using  which  the  plunge  motion  is  then  constructed  with  eqns.  17  and  18.  Recalling  that  LEV 
initiation  in  the  baseline  case  occurs  at  t*  — 1\  =  0.95,  plunge-combined  kinematics  are  constructed 
such  that  LEV  initiation  is  shifted  to  (i)  t*  —  t\  =  0.5,  and  (ii)  t*  —  t\  =  1.5.  The  values  of 
(h/c)  as  determined  from  the  newton  iteration  for  these  two  cases  (hereafter  called  “design  cases”) 
are  —0.5098  and  0.1933  respectively.  As  expected,  the  LESP  theory  predicts  a  negative  plunge 
(equivalent  to  positive  pitch  angle)  to  advance  LEV  formation  and  a  positive  plunge  to  delay  LEV 
formation.  Figure  5.22  shows  the  baseline  case  (pure  pitch),  and  the  two  design  cases  which  are 
constructed  by  combining  the  baseline  pitch  with  the  two  plunge  motions  shown.  Figure  5.23 
illustrated  the  LESP  histories  for  the  baseline  motion  and  two  design  motions.  As  required,  the 
LESP  values  for  the  design  cases  cross  the  critical  value  at  t*  —  t\  =  0.5  and  t*  —  t\  =  1.5. 

CFD  simulations  were  performed  for  the  two  design  cases  to  validate  the  viability  of  the  LESP 
concept  for  advancing  or  delaying  LEV  formation  according  to  specification.  Table  5.3  shows 
vorticity  plots  for  the  baseline  case  and  the  two  design  cases,  at  several  equally  spaced  intervals 
through  the  respective  motions.  The  plots  confirm  that  LEV  formation  can  indeed  be  advanced 
(as  in  designl  to  t*  =  0.5)  or  delayed  (as  in  design2  in  to  t*  =  1.5)  by  suitably  by  superimposing 
a  suitable  plunge  motion.  Using  the  LESP  concept  and  calculating  suitable  superpositions,  LEV 
formation  for  any  arbitrary  motion  may  be  triggered  at  will  or  suppressed  entirely. 


5.4  Interim  conclusions 

In  this  paper,  an  inviscid  theoretical  method  which  handles  large  amplitudes  and  non-planar  wakes 
is  presented,  and  used  to  derive  the  Leading  Edge  Suction  Parameter  (LESP)  which  is  a  measure 
of  the  suction  at  the  airfoil.  Parametric  studies  with  experiments  and  CFD  are  used  to  rigorously 
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Table  5.3:  Vorticity  plots  from  CFD  at  various  instants  during  the  motion 


Baseline 


Designl 


Design2 


test  the  LESP  hypothesis,  that  there  is  a  critical  value  of  the  LESP  for  a  given  airfoil  and  Reynolds 
number  at  which  LEV  formation  is  initiated. 
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In  conclusion,  the  computational  method  used  was  seen  to  reproduce  the  experimental  loads  and 
flow  features  well.  It  was  seen  that  there  is  a  critical  value  of  the  LESP  for  a  given  airfoil  and 
Reynolds  number  at  which  LEV  formation  is  initiated,  except  for  slow-rate  motions  with  significant 
trailing-edge  flow  separation.  The  value  of  critical  LESP  is  completely  independent  of  kinematic 
parameters  such  as  amplitude  and  rate  of  motion,  pivot  location  and  type  of  motion  (pitch  /  plunge). 
By  pre-determining  the  critical  LESP  for  a  given  airfoil  and  Reynolds  number  it  is  possible  to  predict 
whether  LEVs  will  be  formed  for  any  given  motion  kinematics.  Further,  the  LESP  may  be  used 
in  a  design  approach  to  generate  motion  kinematics  which  would  either  prevent  LEV  formation 
or  generate  LEVs  as  per  aerodynamic  requirements.  The  results  demonstrate  that  the  LESP  is  a 
fundamentally  important,  albeit  simple,  theoretical  parameter  that  governs  LEV  formations. 
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Chapter  6 


Model  Reduction  in  Discrete- Vortex 
Methods  for  Unsteady  Airfoil  Flows 


In  this  chapter,  we  propose  a  method  for  model  reduction  in  discrete-vortex  methods  for  two- 
dimensional  flows  on  airfoils.  Discrete-vortex  methods  have  been  successfully  employed  to  model 
separated  and  unsteady  airfoil  flows.2  Earlier  research2  revealed  that  a  parameter  called  the  Lead¬ 
ing  Edge  Suction  Parameter  (LESP)  can  be  used  to  modulate  leading-edge  vortex  (LEV)  shedding 
in  unsteady  flows.  The  LESP  is  a  measure  of  suction  developed  at  the  leading  edge,  and  whenever 
the  LESP  exceeds  a  critical  value,  a  discrete  vortex  is  released  from  the  leading  edge  so  as  to  keep 
the  LESP  at  the  critical  value.  Although  the  method  was  successful  in  predicting  the  forces  on 
and  the  flow  field  around  an  airfoil  in  unsteady  vortex- dominated  flows,  it  was  necessary  to  track 
a  large  number  of  discrete  vortices  in  order  to  obtain  the  solution.  The  current  study  focuses  on 
obtaining  a  model  with  a  reduced  number  of  discrete  vortices  to  model  LEVs,  thus  improving  the 
computation  time.  Vortex  shedding  from  the  leading  edge  is  modeled  by  a  shear  layer  that  com¬ 
prises  a  few  discrete  vortices,  and  a  single  concentrated  vortex  whose  strength  varies  with  time. 
The  single  vortex  at  the  end  of  the  shear  layer  accounts  for  the  concentrated  vortical  structure  that 
comprises  several  discrete  vortex  elements  in  conventional  vortex  methods.  A  merging  algorithm 
is  initiated  when  the  edge  of  the  shear  layer  starts  rolling  up.  Suitable  discrete  vortices  are  iden¬ 
tified  using  a  kinematic  criterion,  and  are  merged  to  the  growing  vortex  at  every  time  step.  The 
reduced  order  method  is  seen  to  bring  down  the  number  of  discrete  vortices  required  to  model  LEV 
structures  significantly.  The  effort  described  in  this  chapter,  although  intended  for  two-dimensional 
flows,  provided  important  algorithms  for  the  UVLM-based  prediction  methods  for  finite- wing  LEV 
formation,  described  in  Chapter  9. 


6.1  Background:  The  LESP-modulated  Discrete  Vortex  Method 
(LDVM) 

This  section  briefly  describes  the  2D  low-order  LDVM  that  uses  discrete-vortex  shedding  from 
the  airfoil  leading-edge  to  simulate  the  intermittent  shedding  of  vorticity  from  an  airfoil  leading 
edge.  The  method  is  developed  by  augmenting  the  large-angle  thin-airfoil  theory,  described  in 
Section  5.1.1.  More  details  of  the  LDVM  approach,  results,  and  comparison  with  computational 
and  experimental  results  are  shown  in  Ref.2 
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6.1.1  LESP  criterion  for  initiation,  shedding,  and  termination  of  LEV 

It  has  been  shown  by  Ramesh  et.  al2  that  the  initiation  and  termination  of  leading-edge  vorticity 
shedding  is  determined  by  a  nondimensional  measure  of  the  suction  at  the  leading  edge.  It  has  also 
been  shown  that  this  leading-edge  section  parameter  (LESP)  is  the  same  as  leading  coefficient  Aq 
in  the  Fourier  series  representing  the  bound  vorticity.  Thus, 

LESP(t )  =  A0(t)  (1) 

With  the  LESP  criterion,  a  given  airfoil  has  a  critical  value  of  LESP  for  a  given  Reynolds  number. 
This  value  is  independent  of  motion  kinematic  parameters  such  as  amplitude,  reduced  frequency, 
and  pivot  location.  Once  the  critical  value  of  LESP  of  a  given  airfoil  is  obtained  for  a  given  Reynolds 
number  using  data  from  experiment  or  CFD  for  a  sample  motion,  LDVM  can  predict  the  initiation 
and  termination  of  intermittent  leading-edge  vorticity  in  any  arbitrary  motion  of  the  airfoil  at  that 
Reynolds  number. 

The  instant  at  which  the  magnitude  of  LESP  exceeds  the  critical  value  marks  the  initiation  of 
leading  edge  vorticity  shedding.  A  discrete  leading  edge  vortex  (LEV)  is  shed  at  every  time  step 
after  initiation  of  leading  edge  vorticity  shedding,  the  strength  of  the  vortex  at  a  time  step  is 
determined  so  as  to  keep  the  LESP  at  its  critical  value.  If  the  LESP  is  positive,  the  sign  of  the 
discrete  vortex  shed  at  the  leading  edge  is  clockwise,  and  the  discrete  vortex  tends  to  get  convected 
towards  the  upper  surface.  On  the  other  hand,  when  the  LESP  is  negative, the  sign  of  the  discrete 
vortex  is  counter  clockwise,  and  it  gets  convected  towards  the  lower  surface.  Vorticity  shedding 
terminates  when  the  magnitude  of  LESP  comes  below  the  critical  value. 

We  follow  the  vortex  placement  method  followed  by  Ansari  et  al.51  to  determine  the  position  of 
the  discrete  vortex  (DV)  generated  at  either  edge  at  any  instant  of  time.  It  is  placed  at  one  third 
of  the  distance  between  the  shedding  edge  and  the  previously  shed  vortex.  The  first  DV  shed  when 
LEV  formation  starts  is  placed  using  the  velocity  at  the  leading  edge.  Like  other  discrete  vortex 
methods,  the  point  vortices  are  convected  with  the  flow  velocity  at  the  location  of  the  vortices.  A 
first-order  time  stepping  process  is  used  since  the  change  in  accuracy  was  not  much  when  high-order 
methods  were  used.  Using  a  vortex  model  with  core  radius  vcore ,  we  have  from  Vatistas  et  al.85 
that  the  induced  velocities  ( u  and  w)  by  the  kth  vortex  in  the  X  and  Z  direction  are: 


Details 

al.2 


2tt  y((X  -  X,)2  +  (Z  -  Zk) 2)2  +  ^ 

in  _  _7 k _ x  -  Xk _  ^ 

2?r  \/((.Y  *  A-,)'-'  i  (Z  Zkff  •  , 

of  force  and  moment  calculations  in  the  discrete  vortex  method  can  be  found  in  Ramesh  et 


6.2  The  reduced  order  model 

This  reduced  order  model  is  developed  to  phenomenologically  model  the  evolution  of  leading  vortex 
structures  as  revealed  by  experimental  and  high  fidelity  CFD  studies.  In  the  early  stages  of  leading- 
edge  vorticity  shedding,  a  shear  layer  is  ejected  from  the  leading  edge.  The  tip  of  this  shear  layer 
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rolls  up  into  a  core  after  some  time.  Once  the  roll-up  is  initiated,  the  shear  layer  keeps  feeding 
vorticity  into  the  core  until  either  the  vorticity  shedding  stops,  or  the  core  pinches  off  from  the 
shear  layer. 

LDVM  predicts  this  flow  pattern  using  discrete  vortices  shed  from  the  leading  edge  for  the  period 
of  time  during  which  the  LESP  exceeds  the  critical  value. The  objective  of  the  reduced  order  model 
is  to  reduce  the  discrete  vortex  count  in  LDVM  by  merging  suitable  pairs  of  discrete  vortices  shed 
from  the  leading  edge.  In  particular,  the  large  number  of  discrete  vortices  that  represent  a  huge 
concentrated  vortex  core  can  be  replaced  by  one  single  vortex  that  has  the  combined  strength  of 
all  the  discrete  vortices. 

We  are  interested  in  studying  the  growth  of  the  core  through  a  model  that  incorporates  a  shear 
layer,  and  a  concentrated  vortex  at  its  tip.  There  are  two  steps  in  modeling  this  phenomenon. 
The  first  is  the  identification  of  the  rollup  of  the  shear  layer.  Once  the  shear  layer  starts  rolling 
up,  a  vortex  that  grows  in  strength  with  time  is  introduced  at  the  tip  of  the  shear  layer.  At  each 
time  step,  some  vorticity  is  added  to  this  growing  single  LEV  (referred  to  as  SLEV  henceforth)  by 
merging  an  appropriate  DV  to  it.  This  procedure  results  in  a  few  point  vortices  forming  a  shear 
layer,  from  which  vorticity  is  fed  into  the  growing  SLEV. 

6.2.1  Initial  rollup 

It  is  proposed  that  the  shear  layer  starts  rolling  up  when  there  is  a  relative  rotation  between  any 
two  successive  vortices  near  the  tip  of  the  shear  layer.  For  identifying  the  onset  of  rollup,  we  keep 
track  of  the  angular  velocity  of  successive  vortex  pairs  in  a  shear  layer.  When  the  angular  velocity 
between  the  ith  and  (i  —  l)th  vortices  in  the  shear  layer  exceeds  a  threshold  value,  the  shear 
layer  is  assumed  to  have  started  rolling  up  at  its  tip.  The  ith  LEV  is  declared  as  the  SLEV. 


6.2.2  Merging  algortihm 

As  observed  by  researchers  like  Moore,31  Sarpkaya17  and  Cortelezzi,30  when  two  like  vortices  come 
near  each  other,  they  start  rotating  about  the  centroid.  If  one  is  significantly  stronger  than  the 
other,  the  weaker  one  appears  to  be  revolving  about  the  stronger  one  since  the  centroid  would 
nearly  coincide  with  the  location  of  the  stronger  vortex.  In  either  case,  the  separation  between  the 
vortices  keeps  decreasing,  and  tends  to  a  steady  state  value  in  the  absence  of  influence  from  any 
other  external  agents  or  boundaries.  Essentially,  they  behave  as  a  system  that  can  be  approximated 
by  a  single  vortex  at  the  centroid,  with  a  strength  that  is  equal  to  the  combined  strength  of  the 
two  vortices.  Following  these  observations,  and  deriving  insight  from  the  flow  pattern  predicted 
by  LDVM,  we  propose  a  merging  algorithm  that  uses  the  angular  velocity  of  the  line  joining  the 
SLEV  and  a  discrete  vortex  in  the  LEV,  the  component  of  relative  velocity  along  this  line  and  the 
separation  between  them  as  the  metrics  for  determining  if  the  discrete  LEV  has  to  be  merged  with 
the  SLEV. 

The  relative  velocity  of  two  vortices  along  the  line  joining  can  be  obtained  in  terms  of  the  relative 
distance  vector  R21  connecting  their  positions,  and  their  relative  velocity  V21  as, 


Vrei  —  V21-R21 

The  angular  velocity,  CI21  between  two  point  vortices  can  be  expressed  as, 


(4) 


45 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


(5) 


^21 


Vn.Rk 

l^il 


where  Rfi  is  the  vector  normal  to  ifei. 

At  every  time  step  after  rollup  is  identified,  the  algorithm  finds  a  DV  within  a  specified  distance 
| /fell  from  the  SLEV  that  has  the  maximum  angular  velocity  Ll  about  the  SLEV  or  the  maximum 
approach  speed  Vrei  towards  the  SLEV.  This  DV  is  merged  to  the  SLEV,  and  the  combined  discrete 
vortex  is  convected  with  the  velocity  induced  at  its  new  location. 

The  search  is  necessary  in  the  initial  phase  of  rollup  when  there  can  be  a  quite  a  few  discrete 
vortices  that  start  rotating  about  the  SLEV.  At  this  stage,  the  shear  layer  may  even  start  rolling 
up  very  close  to  the  leading  edge.  However,  it  has  been  observed  that,  within  a  few  time  steps  after 
rollup  is  initiated,  the  straight  shear  layer  and  the  concentrated  vortex  structure  become  clearly 
distinguishable.  At  this  stage,  the  search  for  a  suitable  merge  candidate  is  not  necessary  anymore. 
Vorticity  can  be  fed  into  the  SLEV  by  merging  the  DV  at  the  tip  of  the  shear  layer. 

We  have  used  a  model  with  a  common  vortex-core  radius  for  all  the  discrete  vortices,  to  prevent  the 
induced  velocities  from  becoming  unbounded  when  two  point  vortices  come  close  to  each  other.2 
The  distance  used  for  searching  in  the  merging  algorithm  has  been  set  by  trial  and  error  to  a  value 
that  is  10  times  this  core  radius.  The  threshold  angular  velocity  to  identify  initial  rollup  was  set  to 
0.001  rad/s.  Also,  it  was  observed  that  the  search  in  the  merging  algorithm  could  be  stopped  after 
10  discrete  vortices  were  merged  to  the  SLEV. 

After  this,  the  DV  at  the  tip  of  the  shear  layer  is  merged  to  the  SLEV  at  every  time  step.  The 
number  of  DVs  in  the  shear  layer  is  approximated  by  the  relation: 


'H'shear 


distance  of  SLEV  from  LE 
0.75  *  core  radius  of  the  DVs 


(6) 


If  the  number  of  DV’s  in  the  shear  layer  is  greater  than  this  value  at  any  time  step,  the  DV  at  the 
tip  of  the  shear  layer  is  merged  to  the  SLEV.  This  replicates  the  shear  layer  feeding  vorticity  into 
the  vortex  core. 


6.2.3  Location  of  the  combined  DV 

When  two  DVs  of  strengths  T 1  and  T2  located  at  X\  and  X2  respectively  are  merged,  it  has  been 
customary  to  add  their  strengths  and  place  the  combined  vortex  at  the  centroid  of  the  two  DVs 
given  by: 


A centroid 


r  1X1  +  r2x2 
Ti  +  r2 


(7) 


However,  this  approach  induces  errors  in  the  flow  field  especially  if  the  merging  takes  place  near 
the  airfoil.30  Cortelezzi  (1992)  proposed  a  merging  scheme  in  which  the  combined  vortex  is  placed 
at  a  different  location  from  the  centroid  so  as  to  conserve  the  velocity  at  the  shedding  location. 
Following  this  idea,  we  propose  a  novel  merging  scheme  in  which  the  merged  DV  is  placed  at  a 
location  referred  to  as  the  optimal  location,  (; xopt,zopt ).  By  placing  the  combined  vortex  at  this 
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location,  we  ensure  that  the  Fourier  coefficients  A$  and  A\  are  conserved  before  and  after  merging. 
This  conserves  the  leading-edge  suction,  as  well  as  the  bound  circulation  of  the  airfoil  given  by 


Tb  =  7 tcU 


Aq  + 


Ai 

2 


A  2D  Newton  iteration  is  employed  to  determine  (xopt,  zopt)  that  conserves  Aq  and  A\. 


(8) 


6.2.4  Identifying  pinch-off 

In  some  cases,  the  shear  layer  stops  feeding  vorticity  into  the  LEV  when  it  attains  a  maximum 
strength.  This  is  followed  by  the  LEV  getting  convected  downstream,  and  the  shear  layer  starting 
to  rollup  near  the  LE  to  form  a  new  LE  vortex  structure.  We  refer  to  this  as  LEV  ’pinch-off’. 
To  model  this  phenomenon,  we  monitor  the  DVs  in  the  shear  layer  at  every  time  step.  If,  at  any 
time  step,  two  adjacent  DVs  in  the  shear  layer  have  opposite  signs  for  their  angular  velocities  with 
respect  to  the  SLEV,  a  pinch-off  is  declared  to  be  identified.  Following  this,  the  DV  in  this  vortex 
pair  that  is  close  to  the  LE  is  declared  as  a  new  SLEV.  The  search  algorithm  is  again  initiated  to 
find  matching  DVs  to  merge  with  this  SLEV.  Meanwhile,  the  DVs  at  the  edge  of  the  shear  layer 
associated  with  the  first  SLEV  are  merged  to  it,  one  at  every  time  step.  Thus,  the  first  SLEV 
represents  the  LEV  structure  that  pinches  off  from  the  shear  layer,  and  the  new  SLEV  models  the 
vortex  core  that  starts  rolling  up  near  the  LE. 


6.2.5  Intermittent  LEV  shedding 

A  new  LEV  is  shed  each  time  the  LESP  reaches  the  critical  value.  This  is  referred  to  as  intermittent 
LEV  shedding.  The  reduced  order  model  handles  this  by  initiating  a  new  SLEV  whenever  LESP 
reaches  its  critical  value.  Each  SLEV  and  its  shear  layer  capture  the  concentrated  core  and  the 
shear  layer  structure  associated  with  the  individual  LEVs. 


6.3  Results 

In  this  section,  we  present  a  comparison  of  the  results  of  the  reduced  order  model  with  those  of  the 
LDVM  for  three  cases. 

6.3.1  Case  1:  Flat  plate  undergoing  0-90  pitching  motion  about  LE 

In  Case  1,  we  study  a  flat  plate  that  is  pitching  about  its  leading  edge  at  a  Reynolds  number  of 
1,000.  The  smoothed  ramp  motion  is  generated  using  Eldredge’s  canonical  formulation.86  At  the 
end  of  the  motion,  the  pitch  angle  reaches  a  value  of  90  degrees.  A  non  dimensional  pitch  rate 
value  of  K  =  0.2  is  used  for  this  motion.  The  critical  value  of  LESP  for  a  flat  plate  at  a  Re  of  1,000 
is  0.11. 2  Fig.  6.1  shows  the  comparison  of  LESP  variation  predicted  by  the  two  methods,  along 
with  the  coefficients  of  lift,  drag  and  moment.  A  comparison  of  the  streamline  patterns  predicted 
by  the  two  methods  at  five  different  instants  of  time  are  shown  in  6.2. 

In  this  case,  the  LESP  reaches  its  positive  critical  value  in  LDVM  during  the  pitch-up  phase  of 
motion,  as  can  be  seen  from  Fig.  6.1.  This  marks  the  onset  of  leading-edge  vortex  shedding  from 
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0.1 


Figure  6.1:  Case  1:  Comparison  of  predictions  of  the  reduced  order  method  and  LDVM.  Variation 
with  t*  =  tU/c  of:  (a)  LESP,  (b)  lift  coefficient,  (c)  drag  coefficient,  and  (d)  pitching  moment 
coefficient  about  the  mid  chord. 


the  upper  surface  of  the  plate.  The  LESP  remains  at  the  critical  value  till  the  end  of  the  motion. 
This  results  in  an  uninterrupted  shedding  of  leading  edge  vorticity  during  the  remaining  span  of 
motion.  The  evolution  of  this  LEV  structure  in  the  LDVM  can  be  observed  in  the  streamline  plots 
in  Fig.  6.2.  The  DVs  representing  the  LEV  are  also  marked  in  these  plots  using  red  dots.  The 
shear  layer  ejected  from  the  leading  edge  soon  rolls  up  into  a  concentrated  vortical  structure  on 
the  upper  surface  near  the  leading  edge.  This  structure  continues  to  grow  in  size  as  the  shear  layer 
keeps  rolling  up  and  feeding  vorticity  into  it.  A  thin  shear  layer  that  connects  the  concentrated 
core  to  the  leading  edge  is  apparent  soon  after  the  shedding  starts.  This  shear  layer  elongates  when 
the  concentrated  vortex  core  convects  away  from  the  plate.  The  concentrated  LEV  core  is  attached 
to  the  LE  through  this  thin  shear  layer  until  the  end  of  the  motion. 

The  reduced  order  model  exactly  predicts  the  onset  of  LEV  shedding  and  thereafter  closely  repli¬ 
cates  the  flow  field.  This  can  be  seen  from  the  streamline  patterns  of  the  reduced  order  method 
in  Fig.  6.2.  An  SLEV  is  initiated  as  soon  as  the  relative  rotation  between  a  pair  of  DVs  at  the 
tip  of  the  shear  layer  exceeds  the  threshold  value.  Following  this,  one  DV  is  merged  to  the  SLEV 
at  every  time  step  using  the  search  criterion.  After  10  DVs  are  merged  to  the  SLEV,  the  search 
is  terminated,  and  the  merging  is  done  based  on  the  approximate  shear  layer  length  given  by  (6). 
By  this  time,  the  thin  shear  layer  and  the  concentrated  core  structure  of  the  LEV  is  apparent  in 
LDVM  (t*  =  1.3).  During  the  later  time  steps,  the  SLEV  in  the  reduced  order  model  represents  the 
strength  and  position  of  concentrated  vortex  core  of  LDVM.  The  elongated  shear  layer  connects 
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the  SLEV  to  the  leading  edge. 

The  comparison  of  the  results  show  that  lift  and  drag  coefficient  histories  as  well  as  the  history  of 
coefficient  of  moment  about  the  mid-chord  predicted  by  the  reduced  order  model  agree  well  with 
those  predicted  by  LDVM.  At  the  last  time  step,  there  are  380  DVs  in  LDVM  that  accounts  for 
the  vorticity  shed  from  the  leading  edge  during  the  period  of  time  for  which  LESP  was  above  its 
critical  value.  Meanwhile,  the  reduced  order  model  can  predict  this  flow  field  by  using  just  89  DVs. 
Corresponding  to  the  decrease  in  the  number  of  discrete  LEVs  in  the  flowfield,  the  run  time  for  the 
reduced  order  (20.1  s)  is  smaller  compared  to  that  of  LDVM  (29.3  s). 

6.3.2  Case  2:  SD7003  airfoil  undergoing  0-90-0  pitching  motion  about  trailing 
edge 

In  this  study,  we  consider  an  SD7003  airfoil  undergoing  a  pitch  up-return  motion.  It  undergoes  a 
smoothed  0-90-0  pitch  up-return  motion  with  the  pivot  point  at  the  trailing  edge.  A  non  dimensional 
pitch  rate  of  K  =  0.4  is  used  for  this  motion.  The  critical  value  of  LESP  for  an  SD7003  airfoil  at 
a  Reynolds  number  of  100,000  is  0.14.  The  comparison  of  LESP  and  the  time  variations  of  force 
and  moment  coefficients  between  the  two  methods  are  shown  in  Fig.  6.3.  Figure  6.4  presents  a 
comparison  of  the  streamline  patterns  predicted  by  the  two  methods  at  different  instants  of  time. 

The  LESP  value  for  this  motion  reaches  the  critical  value  during  the  pitch-up  phase,  as  shown  in 
Fig.  6.3.  It  remains  at  the  critical  value  for  some  time,  and  comes  down  during  the  return  phase. 
This  indicates  that  vorticity  is  shed  from  the  leading  edge  from  the  middle  of  the  pitch-up  phase 
to  the  middle  of  the  return  phase.  The  LESP  value  reaches  the  negative  critical  value  towards 
the  end  of  the  motion,  indicating  the  shedding  of  another  LEV  from  the  lower  surface.  Thus,  this 
is  a  case  where  intermittent  LEV  shedding  occurs.  Fig.  6.4  shows  the  evolution  of  the  flow  field 
predicted  by  LDVM  for  this  case.  The  shear  layer  associated  with  the  first  LEV  starts  rolling  up 
near  the  LE.  This  vortex  continues  to  grow  until  the  vorticity  shedding  stops.  When  the  shedding 
terminates,  the  shear  layer  detaches  from  the  leading  edge.  The  vortex  structure  and  the  shear 
layer  get  convected  downstream,  and  they  interact  with  the  trailing  edge  vorticity.  Meanwhile, 
vorticity  shedding  starts  again  at  the  leading  edge  on  lower  surface  after  a  short  pause.  This  shear 
layer  can  be  seen  to  be  rolling  up  near  the  leading  edge  towards  the  end  of  the  motion. 

Comparison  of  LESP  histories  of  the  two  methods  show  that  the  reduced  order  model  predicts  the 
onset  and  termination  of  leading  edge  vorticity  shedding  at  the  same  instants  of  time  as  LDVM. 
The  reduced  order  model  captures  the  above  mentioned  features  of  the  flowfield  as  can  be  seen 
from  the  streamline  plots  in  Fig.  6.4.  It  predicts  the  roll-up  of  the  first  shear  layer  successfully,  and 
models  the  subsequent  growth  of  the  leading  edge  vortex  very  closely  to  the  predictions  of  LDVM 
using  an  SLEV.  The  SLEV  remains  attached  to  the  leading  edge  through  a  shear  layer  of  varying 
length  airfoil  until  the  vortex  shedding  stops.  During  this  period  of  time,  the  shear  layer  is  seen 
to  add  vorticity  to  the  SLEV.  Upon  the  termination  of  vorticity  shedding,  the  shear  layer  detaches 
from  the  leading  edge,  and  the  DVs  in  the  shear  layer  continue  to  be  merged  with  the  SLEV  in  the 
following  time  steps.  When  the  LESP  reaches  the  critical  value  for  the  second  time,  a  new  SLEV 
is  initiated  to  capture  the  new  LEV  structure.  Meanwhile,  the  DVs  in  the  shear  layer  associated 
with  the  first  SLEV  are  being  merged  with  the  first  SLEV. 

The  force  coefficient  histories  of  the  two  methods  are  in  good  agreement,  except  for  the  slight 
discrepancy  towards  the  end  of  the  motion  when  the  interaction  between  LEV  and  trailing  edge 
vorticity  takes  place.  The  coefficient  of  moment  about  the  quarter  chord  predicted  by  the  reduced 
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Figure  6.2:  Case  1:  Streamline  patterns  predicted  by  LDVM  (left)  and  the  reduced  order  method 
(right)  at  five  time  instants. 
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(a) 


(c)  (d) 


Figure  6.3:  Case  2:  Comparison  of  predictions  of  the  reduced  order  method  and  LDVM.  Variation 
with  t*  =  tU/c  of:  (a)  LESP,  (b)  lift  coefficient,  (c)  drag  coefficient,  and  (d)  pitching  moment 
coefficient  about  the  quarter  chord. 


order  method,  however,  shows  slight  discrepancy  compared  to  the  prediction  of  LDVM.  At  the  end 
of  the  motion,  the  reduced  order  model  predicts  the  leading  edge  vorticity  in  the  flow  field  using 
just  36  DVs  whereas  the  LEV  structure  in  LDVM  is  represented  by  274  DVs.  Computation  time 
savings  in  using  the  reduced  order  model  is  evident  in  the  run  time  of  6.9  s  for  the  reduced  order 
model  as  compared  to  10.6  s  for  LDVM. 

6.3.3  Case  3:  Flat  plate  undergoing  0-45  pitch  up-hold  motion  about  leading 
edge 

This  case  study  deals  with  the  pitch  up-hold  maneuver  of  a  flat  plate  at  a  Reynolds  number  of 
1,000.  The  pitch-ramp  motion  is  generated  using  an  Eldredge  function.  The  flat  plate  pitches  up 
about  the  leading  edge  with  a  non-dimensional  pitch  rate  of  K  =  0.4.  The  value  of  critical  LESP  for 
this  case  is  0.11.  The  comparison  of  LESP  and  the  time  variation  of  force  and  moment  coefficients 
between  the  two  methods  are  shown  in  Fig.  6.5.  Streamline  patterns  predicted  by  the  two  methods 
at  different  instants  of  time  are  compared  in  Fig.  6.6. 

The  LESP  for  this  motion  reaches  the  critical  value  when  the  pitch  angle  is  approximately  3  degrees. 
From  this  instant,  it  stays  at  the  critical  value  till  the  end  of  the  motion.  Consequently,  vorticity  is 
continuously  shed  from  the  leading  edge  from  this  instant  till  the  end  of  the  motion.  The  flow  field 
predicted  by  LDVM  is  shown  in  Fig.  6.6.  The  shear  layer  starts  rolling  up  near  the  leading  edge, 
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25 


Figure  6.4:  Case  2:  Streamline  patterns  predicted  by  LDVM  (left)  and  the  reduced  order  method 
(right)  at  five  time  instants. 
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Figure  6.5:  Case  3:  Comparison  of  predictions  of  the  reduced  order  method  and  LDVM.  Variation 
with  t*  =  tU/c  of:  (a)  LESP,  (b)  lift  coefficient,  (c)  drag  coefficient,  and  (d)  pitching  moment 
coefficient  about  the  quarter  chord. 


forming  a  cluster  of  discrete  vortices.  As  time  progresses,  the  shear  layer  feeds  more  DVs  into  this 
core,  and  the  core  grows  in  size.  The  core  is  seen  to  move  away  from  the  surface  of  the  flat  plate 
towards  the  second  half  of  the  motion.  Later,  the  concentrated  LEV  structure  pinches  off  from  the 
shear  layer.  A  part  of  the  shear  layer  gets  convected  downstream  with  the  core,  while  at  the  same 
time  rolling  up  into  the  core.  The  remaining  part  of  the  shear  layer  that  is  still  attached  to  the 
leading  edge  starts  rolling  up  and  forms  another  concentrated  vortex  structure  near  the  LE. 

It  can  be  observed  from  Fig.  6.6  that  the  reduced  order  model  replicates  the  initial  rollup  and 
growth  of  the  LEV  as  well  as  its  pinch-off.  An  SLEV  is  initiated  at  the  instant  when  the  initial 
rollup  is  identified.  Point  vortices  from  the  shear  layer  are  merged  to  this  SLEV  at  each  time 
step,  and  hence  the  shear  layer  and  the  concentrated  vortex  structure  are  captured  by  the  reduced 
order  model.  The  reduced  order  model  also  successfully  identifies  the  LEV  pinch-off  and  initiates  a 
second  SLEV  to  model  the  rollup  of  the  tip  of  shear  layer  attached  to  the  leading  edge.  Following 
this, the  DVs  in  the  shear  layer  associated  with  the  first  SLEV  are  merged  to  it  at  every  time  step. 
Simultaneously,  the  ones  in  the  shear  layer  associated  with  the  second  SLEV  are  merged  to  the 
second  SLEV.  However,  the  reduced  order  model  predicts  pinch-off  a  little  earlier  compared  to 
LDVM.  This  can  be  noted  from  the  streamline  plot  at  t*  =  5.1.  At  this  time  instant,  a  pinch-off 
is  not  yet  visible  according  to  the  predictions  of  LDVM;  but  a  pinch-off  is  seen  to  be  taking  place 
according  to  the  reduced  order  model.  A  significant  interaction  of  the  trailing  edge  vorticity  with 
the  LEV  core  precedes  the  pinch-off  in  LDVM.  This  may  need  to  be  taken  care  of  for  eliminating 
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Figure  6.6:  Case  3:  Streamline  patterns  predicted  by  LDVM  (left)  and  the  reduced  order  method 
(right)  at  five  time  instants. 
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Case  No. 

Discrete  LEV  count 

LDVM 

Discrete  LEV  count 

reduced  order 

Run  time 
LDVM  (seconds) 

Run  time 

reduced  order  (seconds) 

1 

380 

89 

29.3 

20.1 

2 

274 

36 

10.6 

6.9 

3 

260 

38 

17.1 

13.6 

Table  6.1:  Comparison  of  run  times  for  the  three  cases  between  LDVM  and  the  reduced  order 
method 


the  discrepancy  in  the  predictions  of  the  reduced  order  model. 

It  can  be  observed  from  the  LESP  and  force  and  moment  coefficient  history  plots  in  Fig.  6.5  that 
the  reduced  order  model  accurately  predicts  these  quantities  till  there  is  significant  interaction 
of  leading  edge  vorticity  with  the  LEV  in  addition  to  replicating  the  flow  field.  However,  there 
are  slight  discrepancies  after  the  interaction  becomes  strong.  This  is  evident  in  the  streamline 
plots.  The  trailing  edge  vorticity  starts  interacting  with  the  LEV  after  t*  =  3.6.  The  lift  and  drag 
coefficient  predictions  of  the  reduced  order  model  show  mismatch  approximately  after  this  instant 
of  time. 

At  the  end  of  the  motion,  LDVM  has  260  DVs  in  the  simulation  for  modeling  the  leading  edge 
vorticity  whereas  the  reduced  order  model  has  just  38  DVs.  Owing  to  the  reduction  in  number 
of  discrete  vortices  compared  to  the  first  two  cases,  a  reduction  in  computational  time  has  been 
obtained  with  the  reduced  order  model  for  this  case.  The  run  time  of  the  reduced  order  model  was 
13.6  seconds  whereas  it  took  17.1  seconds  for  the  LDVM  code  to  generate  the  results. 

Table  6.1  provides  a  comparison  of  performance  of  the  new  reduced  order  model  with  that  of  LDVM 
in  terms  of  the  discrete  vortex  count  and  run  time. 


6.4  Interim  Conclusions 

In  this  chapter,  we  presented  an  approach  to  reduce  the  vortex  count  in  discrete  vortex  methods 
for  modeling  unsteady  aerodynamic  flows.  We  proposed  a  phenomenological  approach  to  modeling 
the  leading  edge  vortex  from  2D  bodies  undergoing  unsteady  motion.  We  proposed  an  algorithm  to 
identify  the  initial  roll-up  of  the  shear  layer  shed  from  the  leading  edge,  and  to  reduce  the  discrete 
vortex  count  by  merging  suitable  pairs  of  discrete  vortices.  The  model  has  the  capability  to  handle 
intermittent  LEV  shedding  and  LEV  pinch-off.  The  reduced  order  model  shows  promise  in  terms 
of  reducing  the  count  of  discrete  vortices  compared  to  the  higher  order  LDVM  model,  and  thus, 
saving  computation  time  significantly,  while  also  giving  predictions  of  forces  and  moments  close  to 
that  of  the  higher  order  method.  The  method  is  found  to  be  more  efficient  in  cases  where  there  are 
large  number  of  LEVs  in  the  higher  order  method. 

The  predictions  of  the  model  deviate  from  that  of  the  higher  order  method  when  there  is  significant 
interaction  between  the  vorticity  shed  from  the  leading  and  trailing  edges.  The  model  needs  to  be 
adapted  to  better  handle  such  situations.  The  algorithm  can  be  implemented  to  reduce  the  discrete 
trailing  edge  vortex  count,  and  achieve  faster  computation  times. 
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Chapter  7 


Unsteady  Vortex  Lattice  Method  with 
Modifications  for  Prediction  of  LEV 
Formation  on  Finite  Wings 


In  this  research,  the  unsteady  vortex  lattice  method  (UVLM)  is  used  as  the  foundation  for  the 
low-order  method.  The  traditional  UVLM,  intended  for  inviscid  analysis  of  finite-wing  systems 
undergoing  unsteady  motions,  is  well  understood  and  is  described  in  text  books.76  In  the  first  part 
of  the  current  work,  described  in  Section  7.1,  the  UVLM  was  adapted  with  minor  extensions  for 
prediction  of  LEV  initiation.  These  extensions  relate  to  (i)  calculation  of  the  spanwise  distribution 
of  the  LESP  and  (ii)  an  optional  separated  tip-flow  model  for  improved  modeling  of  the  wing 
tip  vortex.  In  the  second  part  of  the  current  work,  described  in  Section  7.2,  the  UVLM  has  been 
significantly  modified  to  include  a  spanwise- varying  and  time- varying  vortex  sheet  thatis  shed  from 
the  leading  edge.  The  former  technique  is  referred  to  as  “UVLM”  and  the  later  as  “Modified 
UVLM”  in  the  remainder  of  this  report. 


7.1  UVLM  for  prediction  of  LEV  initiation 

This  section  describes  the  minor  modifications  made  in  the  initial  portion  of  this  research  to  the 
standard  UVLM76  for  prediction  of  LEV  initiation.  Figure  7.1  illustrates  the  discretization  of  an 
example  wing-tail  geometry  into  vortex  rings  for  UVLM  analysis.  For  unsteady  effects,  the  time¬ 
stepping  approach76  models  the  shedding  of  trailing-edge  vorticity  by  adding  a  new  vortex  ring 
(per  strip)  at  each  time  step,  as  illustrated  in  Figure  7.2.  The  following  subsections  discusses  minor 
modifications  to  the  UVLM  to  calculate  the  spanwise  distribution  of  LESP  (Section  7.1.1)  and  the 
addition  of  an  optional  separated  tip-flow  model  (Section  7.1.2). 

7.1.1  LESP  evaluation  for  finite  wing 

Of  particular  interest  in  the  current  work  is  the  determination  of  the  spanwise  variation  of  LESP{y ,  t ) 
along  the  wing  at  every  time  step.  For  this  calculation,  the  strengths  of  the  bound- vortex  filaments 
on  each  chordwise  strip  are  considered.  Instead  of  applying  thin  airfoil  theory  to  UVLM  to  obtain 
LESP ,  previous  work87  derived  an  approximation  of  the  first-order  coefficient  Ag(t); 
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Figure  7.2:  Snapshots  of  TEY  shedding  in  vortex  ring  representation. 
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Ao(t)  — 


i.i3ri(t) 


UooC  [cos-1  (1  -  +  sin  {cos-1  (l  -  }] 

The  Ao(t)  term  for  each  wing  strip  is  taken  as  the  LESP(y ,  t)  for  that  strip  at  the  time  step. 


(1) 


7.1.2  Optional  tip- vortex  model 

Traditional  UVLM  assumes  that  the  tip  vorticity  is  attached  to  the  tip-edge  of  the  wing,  and 
sheds  at  the  trailing  edge  of  the  wing  tip.  At  high  angles  of  incidence,  however,  experimental  and 
computational  flow-visualization  images  (for  example,  Figure  3.5)  show  that  the  tip  vorticity  is  not 
attached  to  the  tip,  but  instead  detaches  and  rolls  up  while  being  convected  in  the  approximate 
direction  of  the  freestream.  While  the  effect  of  the  detached  tip  vortex  does  not  differ  much  from 
that  of  the  attached  tip  vortex  for  most  high- aspect-ratio  wings,  discrepancies  for  low- aspect-ratio 
wings  prompted  the  development  of  a  “separated  tip-flow  model”  in  the  current  implementation  of 
the  UVLM. 


Figure  7.3:  Schematic  description  of  separated  tip  flow  model. 


Figure  7.3  describes  the  scheme  of  tip  flow  model.  A  tip  vortex  ring  is  shed  from  the  side  edge 
of  wing  at  each  time  step.  The  basic  idea  of  separated  tip  flow  model  in  current  UVLM  is  based 
on  a  few  conventional  models  for  slender  wings,  as  typified  by  Levin  and  Katz.88  They  assume 
zero  vorticity  at  the  highly  swept  leading  edge  similar  to  the  Kutta  condition  at  the  trailing  edge. 
Considering  these  analogies,  the  strength  of  tip  vortex  ring  Ttr  can  be  approximated  by  Kelvin’s 
theorem  for  a  spanwise  direction  at  wing  tip,  as  follows. 


r  TR,n  r 


n—  1 

b,jmax 


(2) 


It  is  noteworthy  that  the  analogy  has  some  discrepancies.  For  example,  most  of  slender-wing 
systems  have  sharp  leading  edges  but  the  wing  tip  of  a  typical  wing  is  not  always  sharpened,  which 
does  not  support  the  use  of  a  Kutta-condition-like  situation  at  wing  tip.  However,  modelling  the 
tip  flow  is  beyond  the  scope  of  this  project  and  current  UVLM  employs  this  simple  model  when 
evaluating  the  influence  of  tip  vortices.  For  the  location  of  tip  vortex  shedding,  Levin  and  Katz  also 
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(a)  Attached  Tip  Flow  Model.  (b)  Separated  Tip  Flow  Model. 

Figure  7.4:  Comparison  of  distribution  of  vortex  lattices. 


proposes  to  place  the  newly  shed  vortex  ring  at  the  half  interval  of  translation  distance  UooAt/2 
from  the  tip  edge  and  current  UVLM  adopts  the  same  idea. 

To  avoid  the  significant  difficulties  with  modeling  the  roll  up  of  the  separated  tip  flow  and  likely 
interference  of  the  rolled-up  vortex  sheet  with  the  wing  surface,  the  current  separated  tip  flow 
model  is  based  on  straight  wake  sheet,  as  shown  in  figure  7.4.  Figure  7.4  compares  the  vortex-sheet 
geometries  for  the  attached-tip-flow  model  and  the  separated-tip-flow  model. 


7.2  Modified  UVLM  for  prediction  of  LEV  formation 

7.2.1  Leading  edge  vortex 

Leading  edge  vortex  (LEV)  refers  to  vorticity  shed  from  the  leading  edge.  While  TEV  and  tip 
vortices  are  shed  from  respective  edge  in  every  time  step,  LEV  shedding  in  UVLM  is  implemented 
only  when  LESP(t)  >  LESPcra. 

Modelling  of  LEV  for  a  three-dimensional  unsteady  flow  is  a  challenge  because  there  are  few  prior 
successful  works  using  vortex-ring  representation,  in  contrast  to  the  two-dimensional  LEV  model 
in  which  discrete-vortex  representation  has  large  popularity.  On  the  other  hand,  there  are  several 
successful  implementations  of  LEV  models  for  delta  wings.  In  such  applications,  it  is  well  known 
that  LEV  sheet  distribution  of  slender  wing  results  in  a  pair  of  steady  vortical  structure  from  the 
shedding  edges,  and  these  are  easily  modelled  by  a  pair  of  vortex  cores  representing  the  center  line 
of  the  turn,  as  typified  by  Brown  and  Michael.23  In  another  example,  Kandil  et  al. 89  use  multiple 
vortex  filaments  shed  from  the  leading  edge  of  the  delta  wing  in  place  of  LEV  sheets.  The  other 
example  is  a  hybrid  method  of  discrete- vortex  model  and  vortex  ring  representation.  Katz90  assumes 
a  large  vortex  ring  extending  from  the  leading  edge  to  the  trailing  edge  in  which  the  front  edge  of  the 
ring  behaves  as  a  single  LEV  filament  and  the  aft  edge  performs  as  a  single  TEV  filament.  Although 
these  prior  works  give  a  lot  of  insights,  the  current  UVLM  cannot  directly  employ  these  models 
because  of  their  limited  generality.  Also,  line- vortex  model  will  need  a  new  formulation  to  evaluate 
the  representative  circulation  strength  and  line  distribution  in  three-dimension.  Considering  the 
trade-off  between  the  merit  and  the  difficulties  of  both  models,  vortex-ring  representation  was 
chosen  for  the  LEV  model  in  current  work. 

The  main  issue  of  modelling  LEV  is  the  evaluation  of  the  strength  of  circulation  shed  from  the 
leading  edge.  It  is  an  ambiguous  problem  in  UVLM  because  the  velocity  field  at  the  leading  edge 
becomes  numerically  infinite,  but  it  must  be  a  bounded  value.  The  numeric  singularity  at  the 
leading  edge  results  from  the  postulate  of  inviscid  flow.  Thus,  an  additional  condition  to  determine 
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the  bound  circulation  and  LEV  strength  around  the  leading  edge  is  necessary. 

The  simplest  LEV  model  is  to  assume  the  same  Kutta  condition  at  the  leading  edge.  Imposing 
Kutta  condition  at  the  leading  edge,  the  strength  of  nascent  LEV  element  is  simply  obtained  from 
the  adjacent  wing  bound  circulation  closest  to  the  leading  edge  by  Kelvin’s  theorem  and  LEV  sheet 
is  shed  upstream,  like  TEV  is  shed  downstream.  This  method  is  broadly  used  to  model  LEV  of 
delta  wings  with  sharp  leading  edges,  for  example,  by  Levin  and  Katz.88  Ansari  et  al. ,  in  their 
discrete- vortex  model,  derived  LEV  distribution  from  Kelvin’s  theorem  to  obtain  the  flow  field  which 
enabled  Kutta  condition  (zero  vorticity)  at  the  leading  edge.50  However,  in  inviscid  flow  modelled 
by  vortex-ring  representation,  the  wing  bound  circulation  at  the  leading  edge  cannot  be  defined 
clearly  because  the  leading  edge  is  a  singular  point  and  the  vorticity  at  the  leading  edge  becomes 
infinite,  i.e.  numeric  solution  of  the  leading  edge  circulation  strongly  depends  on  the  resolution 
of  discretization.  In  addition,  application  of  Kutta  condition  to  describe  the  flow  field  around 
the  leading  edge  is  still  an  open  question.  Smith91  pointed  out  that  the  flow  velocity  behind  the 
separation  line  must  be  parallel  to  the  stream  line  because  vortex  sheets  form  the  separated  stream 
surface.  Wu  et  al.  proposed  an  asymptotic  viscous  theory  to  model  boundary  layer  separation.92 
Comprehensive  discussion  about  applicability  of  Kutta  condition  in  various  cases  was  provided 
by  Crighton.93  After  a  comparative  investigation  about  flow  perturbation  relating  to  geometry, 
unsteadiness,  separation,  compressibility,  instability  wave,  and  oncoming  turbulence,  he  concluded 
that  Kutta  condition  would  be  valid  in  a  low  Re  and  St  cases,  but  it  is  still  an  open  matter  about 
whether  to  apply  Kutta  condition  to  the  flow  around  the  leading  edge. 

In  the  current  work,  the  LESP  concept  does  not  explicitly  account  for  a  Kutta-type  condition  in 
LEV  simulation  but  has  a  criterion  to  enforce  the  flow  field  at  the  leading  edge.  The  LESP  concept  is 
intended  for  application  to  rounded- leading-edge  geometries.  In  these  cases,  there  must  be  attached 
flow  turning  around  the  leading  edge  before  LEV  formation,  and  this  flow  generates  the  leading- 
edge  suction.  According  to  the  work  of  Ramesh  et  al .,2  there  should  be  a  finite  vorticity  distribution 
at  the  leading  edge  even  under  LEV  formation  because  it  is  known  that  the  leading-edge  suction 
does  not  vanish  but  converges  to  a  finite  value.  For  example,  according  to  Pinkerton94  the  leading- 
edge  suction  force  in  post-stall  condition,  converted  from  the  raw  pressure  data,  asymptotically 
approached  to  a  finite  value  even  though  the  stagnation  point  retreated.  LESP  is  a  measure  of 
suction  and  it  will  also  give  a  direction  for  the  vorticity  at  the  leading  edge  after  initiation  of 
LEV  formation.  Ramesh  et  al.  established  a  new  LEV  model  that  the  vorticity  field  around  the 
leading  edge  saturates  at  LESP(t )  =  LESPcra ,  and  the  strength  of  shed  LEV  is  determined  by 
the  suction  parameter  LESP(t)  holding  an  appropriate  finite  value  LESPcra.  Considering  the 
permanent  continuity  of  vorticity,  because  the  bound  circulation  can  maintain  the  flow  field  around 
the  leading  edge  until  the  onset  of  LEV  initiation  ( LESP(t )  <  LESPcrit ),  it  is  reasonable  to  think 
that  the  vorticity  field  around  the  leading  edge  preserves  the  flow  ( LESP(t )  =  LESPcrit )  and  the 
surplus  of  vorticity  is  eventually  shed  to  the  fluid  as  LEV.  The  current  UVLM  employs  the  LESP 
concept  including  this  LEV  model.  The  strength  of  shed  LEV  is  determined  using  an  iteration  in 
which  LESP(t)  saturates  at  LESPcrn  for  every  wing  strip  in  which  LESP(t)  exceeds  LESPcra. 
Denoting  superscript  k  as  iteration  number,  the  strength  of  LEV  sheet  for  a  wing  strip  can  be 
evaluated  by  Newton- Raphson  method  as  follows: 

r|  =  LESpllZ^lpk^SPcrit  -  LESPk~l)  +  I*"1  (3) 

Due  to  the  assumption  of  nonzero  vorticity  at  the  leading  edge  under  LEV  formation,  LEV  sheet 
is  tangentially  shed,  which  means  LEV  sheet  leaves  in  the  z  direction  in  body  frame  of  reference. 
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(a)  At  t  =  (n  —  l)At, 
before  LEV  convection. 


^L,n-2 


(b)  At  t  =  nAt, 
before  LEV  calculation. 


^L,n-2 


(c)  At  t  =  nAt, 
after  LEV  calculation, 
before  LEV  convection. 


Figure  7.5:  Placement  of  LEV  ring  vortex. 


In  order  to  place  the  new  LEV  sheet,  the  current  UVLM  employs  a  method  of  LEV  placement 
proposed  by  Ansari  et  al. 51  The  idea  is  to  place  new  point-vortex  of  LEV  at  one-third  distance 
from  the  leading  edge  to  the  last  LEV.  It  can  be  formulated  as  below, 


1 


7^%i,L.E.  +  ^TL,n-2 


(4) 


This  formulation  is  an  approximated  conversion  from  grid-based  vortex  sheet  model  to  discrete 
vortex  model  and  is  advantageous  in  sense  that  it  can  simultaneously  follow  the  wing  motion  and 
vortex  convection. 


As  shown  in  figures  7.5b-7.5c,  the  placement  of  a  new  LEV  panel  is  implemented  by  splitting  a 
previously  shed  LEV  panel  into  two  LEV  panels  at  the  point  determined  by  equation  4.  Then, 
in  figure  7.5,  a  nascent  LEV  T is  defined  as  the  vortex  ring  lying  across  the  leading  edge  and 
%L,n- 1-  However,  as  shown  in  figures  7.5,  the  front  end  of  the  new  LEV  ring  coincides  with  the 
leading  edge  and  this  component  of  the  LEV  ring  violates  the  precedent  Kutta-type  assumption 
because  it  behaves  as  a  large  and  negative  circulation  on  flow  field. 


Leading  edge  vortex  sheet 


In  order  to  prevent  this  numerical  side  effect,  the  current  UVLM  introduces  a  pseudo  vortex  ring  to 
the  distribution  of  the  wing  bound  circulation.  Figure  7.6  shows  the  alignment  and  the  strength  of 
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the  pseudo  vortex  ring.  The  pseudo  vortex  ring  is  a  fictitious  vorticity  distribution  and  compensates 
for  the  mathematical  discontinuity  caused  by  LEV  rings  at  the  leading  edge.  The  pseudo  vortex 
ring  is  overlaid  on  the  wing  bound  circulation  from  separation  line  (the  leading  edge)  to  the  trailing 
edge  shedding  point  (not  always  the  same  as  the  trailing  edge).  Considering  Helmholtz  law  at 
the  leading  edge,  the  strength  of  the  pseudo  vortex  ring  can  be  derived  as  the  current  Tl  at  the 
wing  strip.  The  presence  of  the  pseudo  vortex  ring  affects  the  circulation  conservation  (Kelvin’s 
theorem),  i.e.  it  provides  the  circulation  component  T^  at  the  trailing  edge  shedding  point,  which 
is  accounted  for  in  the  modified  UVLM. 


7.2.2  Solution  of  linear  system 

As  noted  in  previous  section,  the  distribution  of  wing  bound  circulation  is  derived  by  boundary 
condition  on  wing  surface  which  specifies  that  zero  normal  component  of  the  velocity  on  the  surface 
is  enforced.  This  boundary  condition  can  be  described  by  linear  equation  5. 


(vb  +  Vqo  +  vm  +  vw  +  vTl  +  vtr  +  vL)  •  ft  =  0  (5) 

As  mentioned  in  previous  sections,  velocity  induced  by  motion  vm  is  determined  by  given  motion 
kinematics,  velocity  induced  by  shed  vortices  %,  vtl,  ^tr ,  and  vl  are  given  by  application  of 
Biot-Savart  law  to  each  shed  vortex  ring  in  the  flow  field.  Thus,  unknown  values  are  only  wing 
bound  circulation  T&  and  equation  5  is  rewritten  for  K th  wing  bound  vortex  ring  and  the  other 
known  vortices  using  influence  coefficient  as  below. 


N 

^a&,KnT&,n  +  (hoo  +  Vm,K)nK  +  a tiXtl  +  ^^trFtr  +  =  0  (6) 

n= 1 


/rM\ 

^  (Vx>  H”  Vm, l)Rl  ^ 

(^W,  1  ^ 

K,nm] 

^bji 
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(ffoo  +  Vm,K)nK 

Fw,K 

KTb,Nj 

^(Vx)  +  Vm,N)vi]\[  j 

\FW,nJ 
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(  Ttr,i  ^ 

rL)1  \ 

+  [aTL,nm\ 

^TL,K 

+  [ dTR,nm\ 

r  tr,k 

T 

Fl,K 

\Ttr,nJ 

•  •  4 

PH, 

In  equation  7,  the  symbol  [  ]  means  a  N  x  N  matrix.  As  a  result,  the  boundary  condition  for 

wing  bound  circulation  reduces  to  an  N  x  N  matrix  equation,  and  the  unknown  T^  is  derived  by  a 
linear  solution  technique,  such  as  LU  decomposition  method. 
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7.2.3  Calculation  of  aerodynamic  load 

Aerodynamic  loads  are  obtained  from  pressure  distribution  on  wing  surface  and  the  leading-edge 
suction  force.  The  pressure  distribution  is  derived  fron  the  wing  bound  circulation  distribution  and 
the  suction  force  is  evaluated  from  the  Aq  value  for  each  wing  section. 


Cl  =  Cn  cos  a  +  Cs  sin  a 


(8) 


7.2.4  Vortex  sheet  roll  up 

Vortex  roll  up  in  UVLM  is  a  model  for  wake  vorticity  convection.  Vorticity  distribution  moves  by 
the  local  flow  velocity.  Vortex  roll  up  is  important  in  situations  in  which  shed  vortices  directly 
interact  with  the  lifting  body.  Under  LEV  formation,  roll  up  becomes  a  dominant  feature  of 
aerodynamics91  because  the  LEV  structure  develops  near  the  wing  and  generates  vortex  lift.  On 
the  other  hand,  vortex  roll  up  is  often  problematic  due  to  a  numeric  instability.  Characteristic 
numeric  instability  in  UVLM  is  caused  by  singularity  of  vortex  model  and  geometric  discrepancy, 
which  are  discussed  in  the  following  subsections. 

7.2.4. 1  Numerical  process  of  vortex  roll  up 

Vortex  wake  roll  up  calculates  a  displacement  of  distribution  of  vortex  sheet  driven  by  local  velocity 
field  vr.  The  displacement  of  a  grid  point  of  a  vortex  ring  during  a  unit  time-step  in  inertial  frame 
of  reference  Ax  is  simply  described  as 

Ax  =  vrA  t,  (9) 

where  vr  =  (v^ + v\y + + vtr + )  •  Each  velocity  term  is  given  by  Biot-Savart  law,  but  the  finite 
difference  method  (FDM)  for  the  time  advancing  scheme  has  a  number  of  options.  Several  kinds 
of  numerical  schemes  are  applicable  to  UVLM  and  the  simplicity  of  the  scheme  is  proportional  to 
computational  speed  and  numeric  instability.  A  basic  way  to  calculate  the  time-marching  vortex 
roll  up  is  Euler  method  which  estimates  the  displacement  only  by  information  from  the  current 
time-step,76  shown  below. 

Ax  =  v?At  (10) 

The  well-known  multi  time-step  method  to  solve  ODE  is  Runge-Kutta  method.  Runge-Kutta 
method  is  based  on  Taylor  expansion  and  fourth-order  Runge-Kutta  method  is  especially  popular 
in  use; 

Ax=^(yK1  +  K2  +  K3  +  K4) 
where  each  coefficient  vectors  are  available  as 

Ki  =  At  •  vr(t,  x) 
i?2  =  At  •  vr  (t  +  \At,  x  +  \Ki^J 
i?3  =  At  •  vr  (t  +  \  At,  x  + 
i?4  =  At  •  vr  (t  +  \  At,  x  +  K^j  . 


(11) 


(12) 
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However,  Dovgii  and  Shekhovtsov95  pointed  out  that  conventional  multi-step  method  would  likely 
spoil  the  accuracy  in  UVLM.  The  paradox  that  high-order  scheme  enlarges  the  error  results  from 
the  fact  that  major  source  of  error  in  UVLM  is  an  input  error,  and  a  high-order  scheme  with 
multiple  references  increases  the  chance  of  input  error.  For  example,  fourth-order  Runge-Kutta 
method  requires  us  to  calculate  the  derivative  four  times  to  estimate  the  four  unknown  coefficients, 
and  this  process  accumulates  the  input  errors.  In  order  to  prevent  this  paradoxical  error  in  multi- 
step  method,  the  UVLM  has  to  employ  a  method  with  single  reference.  It  is  known  that  some 
multi-step  method  with  variable  coefficients  are  replaced  by  constant  coefficients  by  approximating 
the  reference  function  as  an  algebraic  series.  For  example,  Paul  et  al .96  employed  Adams-Bashforth 
method  in  UVLM. 


In  addition,  according  to  Dovgii  and  Shekhovtsov,95  the  order  of  time-marching  method  in  UVLM 
does  not  determine  the  numeric  accuracy  and  they  proposed  two  low-order  FDMs  as  an  alternate  of 
predictor-corrector  scheme.  One  FDM  that  they  proposed  is  a  two-step  Euler  method  with  second 
order  approximation,  as  shown  below. 


Ax 


=  di  +  ALELri(TkEl 


tn+1  -  V 


(13) 


Equation  13  takes  into  account  an  acceleration  of  vortices  by  putting  weight  in  present  velocity 
term.  Another  FDM  is  a  balanced  FDM  of  equation  13. 


Ax  = 


At 


(14) 


The  FDM  of  equation  14  consists  of  the  present  and  the  previous-time-step  information,  and  it 
is  a  first-order  approximation.  After  investigations  and  tests  of  these  schemes  in  current  UVLM 
problem,  it  was  shown  that  high-order  scheme  gives  less  improvement  in  accuracy.  For  this  reason, 
current  UVLM  takes  the  opinion  of  Dovgii  and  Shekhovtsov95  and  employs  equation  14  as  the 
time-advancing  method  for  vortex  roll  up. 


7. 2.4. 2  Desingularized  vortex  model 


As  seen  in  previous  sections,  UVLM  employs  vortex  ring  (closed  vortex  filament)  discretization. 
However,  a  mathematically-derived  vortex  filament  is  a  singularity  at  the  center  which  could  cause 
extremely  large  velocity  and  physically  impossible  vortex  sheet  intersection.  This  discrepancy 
results  from  the  assumption  of  inviscid  flow  and  vortex  lumping  in  that  all  intensity  of  three- 
dimensionally  distributed  vorticity  is  concentrated  on  a  line.  The  velocity  field  induced  by  vortex 
sheets  obeys  Biot-Savart  law  and  this  has  0(— 1)  order  in  distance  between  vortex  filaments  as 
variable.  When  they  are  close  each  other  and  r  — >  0,  a  great  error  occurs  in  the  magnitude 
of  velocity,  namely  an  input  data  error.95  According  to  an  error  analysis  for  Cauchy  problem 
reported  by  Dovgii  and  Shekhovtsov,95  the  input  errors  in  the  system  of  VLM  become  far  greater 
than  the  truncation  errors.  Thus,  some  viscous  model  is  required  to  eliminate  the  mathematical 
singularity. 


Real  vortex  structure  has  a  finite  vorticity  distribution  around  the  center.97  As  the  vortex  model 
which  guarantees  finite  value  in  vortex  center,  Rankine  model  (equation  15)  and  Lamb-Oseen  model 
(equation  16)  are  well-known. 


Vo  = 


(r  >  rc) 
(r  <  rc) 


(15) 
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(16) 


Both  models  consist  of  the  original  Biot-Savart  law  and  an  additional  nondimensional  cutoff  func¬ 
tion. 


Radius,  m 


Figure  7.7:  Comparison  of  tangential  velocity  distribution  of  Rankine  and  Lamb-Oseen  model 
(Assuming  T  =  1,  rc  =  2 \[vi  =  0.5). 


As  shown  in  figure  7.7,  the  Rankine  vortex  model  has  C 1  discontinuity  at  the  cutoff  radius  rc, 
dividing  internal  Couette  flow  region  and  external  circulation  flow  region.  In  contrast,  Lamb-Oseen 
vortex  model  shows  a  continuous  distribution  with  time-dependent  viscous  diffusion  distance 
(It  should  be  noted  that,  Lamb-Oseen  equation  16  is  known  as  a  particular  solution  of  Navier-Stokes 
equation  for  a  single  vortex  but  not  general  solution  for  multiple  vortex  system.98)  Rankine  vortex 
model  is  0(r3)  approximation  of  Lamb-Oseen  model,  which  is  confirmed  by  the  assumption  that 
rc  =  2 y/vt  and  applying  Taylor  expansion  to  the  cutoff  function  of  equation  15  and  16  as  follows. 


1  me  4 ut  =  \  —  —  -  — 


As  an  example  to  use  Rankine  vortex  model  in  UVLM,  Gennaretti  and  Bernardini"  reported  a 
reasonable  accuracy  in  accordance  with  their  numerical  investigation.  It  is  known  that  the  discon¬ 
tinuity  of  Rankine  model  is  solved  by  adding  a  limit  radius  in  radial  coordinate  r  of  denominator 
of  vortex  model.  Vatistas  et  al. 85  attempted  a  comprehensive  analysis  of  application  of  Rankine 
vortex  model  and  proposed  an  empirical  formula  as  below. 


ve  = 


(17) 


Equation  17  was  utilized  by  Ansari  et  al. 51  and  Ramesh  et  al ?  and  proven  good  performance. 
The  cutoff  function  of  equation  17  is  suitable  for  scalar  form  but  not  always  easy  to  apply  in 
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vector  form.  Thus,  Lamb-Oseen  vortex  model  is  useful  in  current  work.  The  application  of  Lamb- 
Oseen  vortex  model  to  a  vortex  method  was  reported  by  Kuwahara  and  Takami98  and  Summa.100 
According  to  their  study  about  vortex  roll  up,  Lamb-Oseen  vortex  model  performs  good  convergence 
when  v  >  0.002  and  shows  numeric  instability  in  a  smooth  media  v  <  0.002.  This  result  proves 
that  Lamb-Oseen  vortex  model  can  alleviate  the  singularity  of  vortex  model  but  brings  a  new 
arbitrariness  to  UVLM  in  the  sense  that  the  viscosity  of  fluid  v  numerically  influences  the  result 
of  inviscid  flow  model.  About  this  arbitrariness  in  UVLM,  Dovgii  and  Shekhovtsov  proposed 
to  employ  singular  radius  as  the  cutoff  radius  of  the  artificial  viscosity  rc.95  To  contrast  that 
viscous  diffusion  in  Lamb-Oseen  vortex  model  is  evaluated  physically,  singular  radius  is  derived  by 
numerically  expected  precision  determined  by  resolution  of  discretization.  Dovgii  and  Shekhovtsov 
pointed  out  that  it  is  impossible  to  expect  correct  value  from  Biot-Savart  law  at  distance  less 
than  half  of  panel  size,  since  nobody  can  require  a  higher  precision  of  the  model  than  that  given 
by  the  boundary- value  condition  at  the  surface.95  Accordingly,  the  viscous  diffusion  distance  in 
equation  16  can  be  replaced  by  rc  without  numeric  setback.  For  these  reasons,  Lamb-Oseen  model 
restricted  by  singular  radius  is  thought  as  the  most  suitable  cutoff  function  to  desingularize  current 
UVLM  and  Biot-Savart  law  is  modified  as  follows. 


v(r)  = 


r  fsxre 

4tt  \r 3  x  re |2 


m  +  mi) 


(18) 


7. 2. 4. 3  Impingement  of  vortex  sheets 

In  vortex  roll  up,  discrete  vortices  or  grid  of  vortex  rings  tend  to  congest  in  a  place  and  form  spiral 
structure.101  The  cutoff  function  (Equation  18)  mitigates  the  overestimation  of  velocity  induced  by 
proximal  vortex  filament  but  does  neither  inhibit  this  tendency  of  vortex  congestion  and  nor  prevent 
the  eventual  irregularity  in  highly  concentrated  vortex  region.97  Thus,  vortex-sheet  impingement 
could  happen  in  fully  developed  vortex  structure  even  if  the  UVLM  used  a  desingularized  vortex 
model.  For  some  aerodynamic  problem  in  which  the  type  of  wing  is  inevitably  exposed  to  a 
vortical  wake,  such  as  LEV,  rotor  craft,  and  turbine  system,  this  vortex-wing  impingement  matters. 
For  wing- wake  impingement,  several  directions  of  solution  has  been  proposed.  Gennaretti  and 
Bernardini"  solved  the  vortex-blade  interaction  in  rotor  craft  by  UVLM.  In  their  work,  the  trailing 
edge  wakes  were  divided  to  the  near- wake  and  the  far- wake  segment  from  the  shedding  edge. 
According  to  their  investigation,  the  velocity  component  induced  by  the  far-wake  segment  was 
represented  by  blade  bound  circulation  and  the  near- wake  segment.  The  threshold  for  distinguishing 
between  near-wake  and  far-wake,  however,  was  not  clearly  introduced.  Eventually,  their  UVLM 
simulation  showed  the  rotor  blades  immersed  in  vortex  wake  sheets  but  it  does  not  impair  the 
calculation  result.  The  same  problem  was  solved  by  Wie  et  a/.,102  in  which  they  proposed  to 
realign  the  trailing  edge  wake  potential  by  integrating  velocity  field  through  an  arbitrary  path, 
where  the  velocity  field  was  determined  by  original  vortex  distribution.  This  method  does  not 
provide  direct  solution  of  velocity  field  but  performed  reasonably.  These  solutions  are  successful 
but  brings  another  arbitrariness  in  calculation  and  it  is  not  preferable  for  the  present  work.  Then, 
for  current  UVLM,  a  special  treatment  for  unreal  vortex  sheet  penetration  to  wing  is  employed. 
The  idea  of  the  treatment  is  to  assume  a  slip- wall  condition,  which  is  originally  proposed  by  Suresh 
Babu  et  al  (from  a  personal  suggestion).  In  computation  of  vortex  roll  up,  the  displacement  for  a 
grid  of  vortex  sheet  which  is  closer  to  the  wing  surface  than  singularity  radius  rc  was  evaluated  as 
surface  tangential  velocity  as  below. 

Ax  =  (Too  +  Vm  +  VW  +  VTL  +  VTR  +  VL  +  VT)A  t  (19) 
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(a)  Top  view 


(b)  Side  view 

Figure  7.8:  Range  of  application  (colored  by  gray). 


This  solution  forces  vortex  sheet  near  the  wing  surface  to  slide  in  the  tangential  direction  but 
would  interrupt  vortex  shedding  around  the  wing  edges.  Therefore,  the  range  of  application  for  the 
prevention  of  impingement  must  have  a  padding  through  the  edge  lines,  as  shown  in  figure  7.8.  In 
figure  7.8,  gray  regions  denote  the  region  where  equation  19  is  applied,  and  about  20%  of  padding 
in  chord  is  used  in  the  present  work. 


7.2.5  Amalgamation  of  vortex  sheets 


In  the  context  of  vortex  methods,  amalgamation  means  to  merge  some  discrete  vortices.  The 
motivation  of  amalgamation  method  for  vortex  roll  up  are  roughly  categorised  into  physical  reason, 
numeric  reason,  and  economic  reason.  First,  potential  flow  theory  and  UVLM  do  not  simulate 
viscous  dissipation  and  turbulence.  Thus,  the  shed  wakes  are  eternally  preserved  in  the  flow  field. 
In  case  the  preservation  of  vortex  sheets  is  not  preferable,  amalgamation  method  should  be  used. 
Second,  vortex  roll  up  develops  the  vortex  structure  but  it  tends  to  cause  an  irregular  distortion 
and  intersection  in  vortex  sheet  distribution.  Amalgamation  method  can  smoothen  these  geometric 
discrepancies.  Third,  merging  vortex  sheets  reduces  the  number  of  vortex  rings  and  contributes  to 
save  computation  time. 

According  to  Sarpkaya,97  the  idea  of  amalgamation  was  suggested  by  Ham103  to  inhibit  the  error 
in  calculation  of  induced  velocity  by  proximal  discrete  vortices.  Ham  provided  the  basic  scheme  of 
vortex  sheets  amalgamation  and  the  basis  of  the  scheme  has  not  been  largely  changed  by  successor 
methods.  The  main  idea  of  amalgamation  of  vortex  sheets  is  to  give  a  united  circulation  which 
has  the  summation  of  the  strength  of  the  merged  circulations  and  is  located  at  their  centroid. 
Considering  amalgamation  from  Nth  to  Mth  of  discrete  vortices  into  new  one,  the  strength  and 
the  location  of  the  new  representative  discrete  vortex  Tnew  and  x^new  is  formulated  as  follows. 


r 


new 


n=N 


(20) 
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The  formulation  of  amalgamation  method  as  shown  in  equation  20  and  21  must  guarantee  the  far 
field  boundary  flow  condition.104  Although  amalgamation  method  brings  aforementioned  benefits 
in  UVLM,  the  process  is  fundamentally  a  nonlinear  and  irreversible  approximation  in  the  sense 
that,  though  the  circulation  strength  is  conserved,  the  resulting  velocity  field  around  the  region  of 
amalgamated  vortices  is  not  identical  to  that  prior  to  amalgamation.97  Further,  it  is  known  that 
the  merging  vortices  near  the  lifting  body  could  cause  a  discontinuity  in  force  calculation.97  Thus, 
amalgamation  method  should  be  used  with  care. 

Amalgamation  schems  have  been  studied  by  several  engineers  from  various  perspectives.  Spalart104 
provided  a  theoretical  reinforcement  to  the  basic  scheme  as  shown  in  equation  20  and  21  by  Taylor 
expansion  analysis  and  derived  a  merging  criterion  to  minimize  the  truncation  error.  Spalart’s 
criterion  was  largely  simplified  and  employed  in  the  work  of  Ansari  et  al. 51  Moore31  focused  on 
the  geometry  of  the  rolled  up  vortex  structure  and  established  another  amalgamation  criterion 
considering  numeric  resolution.  In  his  investigation,  there  must  be  a  minimum  number  of  discrete 
vortices  to  be  able  to  depict  the  spiral  structure  of  rolled  up  vortex  and  main  numeric  discrepancy 
was  caused  by  low  resolution  to  capture  the  internal  spiral  structure.  Thus,  the  amalgamation 
was  carried  out  in  case  that  the  number  of  grid  of  vortex  rings  cannot  correctly  capture  the  spiral 
distribution.  The  idea  to  solve  the  amalgamation  criteria  from  geometric  viewpoint  was  modified 
by  Lamarre  and  Paraschivoiu,105  in  which  the  highly  rolled  up  portion  was  detected  by  bending 
angle  of  vortex  rings  and  merged.  Lamarre  and  Paraschivoiu  also  proposed  adaptive  panelling 
scheme,  namely  to  compensate  the  lack  of  discrete  vortex  by  reproducing  another  vortex  filament 
in  different  area.  More  intuitive  approach  was  attempted  by  Suresh  Babu  et  al.106  According  to 
the  observation  of  Suresh  Babu  et  al. ,  a  vortex  structure  can  be  represented  by  a  strong  vortex 
core  and  it  is  grown  by  absorbing  the  surrounding  discrete  vortex  components  which  is  detected  by 
relative  angular  velocity. 

For  the  present  work,  the  current  UVLM  implementation  aims  to  mitigate  a  numeric  instability 
and  to  equip  a  geometrically  fundamental  function  to  calculate  the  development  of  free  wake.  The 
numeric  instabilities  in  UVLM  is  typical  due  to  the  intersection.  On  the  one  hand  the  intersection 
of  vortex  and  wing  has  been  introduced  in  previous  sections.  On  the  other  hand  the  object  of 
amalgamation  method  is  to  solve  the  intersection  of  vortex  sheets  itself,  which  is  often  observed 
in  the  spiralling  region  of  rolled  up  vortex  sheets.  The  geometrically  fundamental  function  deform 
the  two-dimensional  vortex  sheets  to  the  three-dimensional  manifold  and  this  problem  is  discussed 
in  the  subsequent  section. 

In  order  to  detect  these  unrealistic  features  in  vortex  sheet  distribution,  current  UVLM  uses  a 
modification  of  past  methods  proposed  by  Moore31  and  Lamarre  and  Paraschivoiu.105  The  ex¬ 
ternal  angle  formed  by  adjacent  vortex  rings  are  evaluated  in  every  time-step  and  if  this  shows 
characteristics  of  adverse  features  like  wave  and  intersection,  the  amalgamation  scheme  is  applied 
to  the  region  at  hand.  The  external  angle  consists  of  three  components:  angle  measured  on  a  XY 
plane,  YZ  plane,  and  XZ  plane.  Denoting  a  pair  of  vectors  a  and  b  on  a  plane,  the  external  angle 
0  which  they  forms  is  calculated  by  an  intrinsic  functions  as  follows. 


6  =  sgn(a  x 


cos 


(22) 


From  discrete  geometry,  the  source  of  vortex  sheets  intersection  can  be  defined  as  pleat  structure 
and  turn  structure.  Turn  structure  forms  spiral  pattern  in  tip  vortex  (assuming  no  secondary 
vortex)  and  wave  pattern  in  TEV  and  LEV  except  for  the  starting  vortex.  Figure  7.9  describes 
these  two  structures.  The  pleat  structure  can  be  detected  as  the  condition  that  the  sign  of  the 
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Figure  7.9:  Schematic  description  of  vortex  sheet  structure  to  be  amalgamated. 


Figure  7.10:  Schematic  description  of  implementation  of  amalgamation. 


external  angle  changes  in  three  times  in  a  row.  The  structural  condition  of  pleat  can  be  described 
as 


On-i  >  0,  On  <  0,  6n+ i  >  0 
6n-i  <  0,  On  >  0,  6n+ i  <  0 


(23) 


Taking  figure  7.9a  as  an  example,  62  <  0,  #3  >  0,  #4  <  0,  the  3rd  grid  of  vortex  sheets  should  be 
merged.  The  turn  inside  of  the  spiral  vortex  sheet  structure  can  be  detected  by  the  region  whose 
summation  of  a  series  of  the  external  angles  exceeds  a  threshold  6a.  If  6a  =  2i r,  the  amalgamation 
scheme  would  merge  the  inner  turn  and  the  only  outermost  circle  would  remain.  It  can  be  described 
by 

M 

n=N 


>0a 


(24) 


and  the  vortex  rings  which  is  from  Nth  to  Mth  should  be  amalgamated. 


The  procedure  of  amalgamation  is  illustrated  in  figure  7.10.  Figure  7.10  shows  that  the  inner  turn 
indicated  by  dashed  line  and  vortex  rings  with  suffix  from  ns  to  ne  + 1  is  merged  into  one  circulation 
denoted  by  suffix  ne  +  1.  The  detailed  procedure  is  explained  in  figure  7.11.  As  mentioned  before, 
UVLM  is  based  on  vortex  ring  discretization  and  physically  conserved  parameter  is  not  the  strength 
of  vortex  ring  Tn  but  the  difference  of  circulation  between  adjacent  vortex  rings  T'n  =  Tn+ 1  —  Tn. 
Therefore,  converting  equation  20  to  the  vortex  ring  discretization,  representative  circulation  at 
the  centroid  must  be  Tns  —  Tne+ 1.  As  shown  in  figure  7.11b,  this  means  that  the  amalgamation 
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N  N-l  N-2  ns+l  ns  ns-l  ns-2  ne+l  nc  ne-l  nc-2  2 


(a)  Detection 
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N  N-l  N-2  ns+l  ns  ne+l  nc  ne-l  ne-2  2  10 

••• 
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(b)  Amalgamation 
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••• 
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(c)  Index  shift 

Figure  7.11:  Schematic  description  of  procedure  of  amalgamation. 
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process  simply  eliminates  the  vortex  rings  from  index  ns  —  1  to  ne  +  2  and  connects  the  nsth  and 
ne  +  1th  vortex  ring  at  the  centroid 


ns  —  1 

5]  xn(rn+1  -  rn) 

n=ne+ 1 
xne+l,new  =  ~  _1 

E  (r„+i  -  r„) 

n=ne+ 1 


ns— 1 

E  £n(rra+i  -  rn) 

n=ne+l 

rns  —  rne+i 


(25) 


After  the  amalgamation  scheme  is  processed,  the  indices  of  vortex  rings  which  are  located  in 
upstream  from  the  amalgamation  centroid  are  shifted  as  shown  in  figure  7.11c. 


7.3  Summary  of  Modified  UVLM 

Figure  7.12  shows  a  flowchart  for  the  modified  UVLM  algorithm.  The  algorithm  for  the  conventional 
UVLM  roughly  consists  of  three  main  steps:  (1)  evaluation  of  the  velocity  induced  by  free  stream, 
motion,  and  wakes,  (2)  derivation  of  the  wing  bound  circulation  distribution  from  surface  boundary 
condition,  and  (3)  calculation  of  the  pressure  distribution  from  the  obtained  velocity  field  on  the 
surface.  The  modified  UVLM  from  the  present  work  includes  new  contributions  for  simulating 
LEV  formation  using  the  LESP  concept  proposed  for  2D  flows  by  Ramesh  et  al. ,2  and  incorporates 
capability  for  flow  simulation  of  LEV-dominated  unsteady  finite  wing  system. 

Results  from  the  UVLM  are  applied  to  prediction  of  LEV  initiation  in  chapter  8  and  to  prediction 
of  LEV  formation  in  chapter  9. 
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Chapter  8 


Low-Order  Prediction  of  LEV 
Initiation  on  Finite  Wings 


This  chapter  presents  the  results  of  studies  on  LEV  initiation  on  finite  wings,  with  emphasis  on 
the  effectiveness  of  the  low-order  prediction  capability  using  the  UVLM  code.  The  studies  in  this 
portion  of  the  effort  were  carried  out  with  the  objective  of  answering  the  following  key  questions: 

(i)  For  any  wing  shape  and  motion,  what  are  the  values  for  time  instant,  angle  of  attack,  and 
spanwise  location  at  which  CFD  predicts  initiation  of  LEV  formation?  How  do  these  values  vary 
with  wing-geometry  and  motion-kinematic  parameters? 

(ii)  For  any  wing  shape  and  motion,  how  does  the  maximum  in  the  UVLM-predicted  spanwise 
variation  in  the  LESP  at  the  time  instant  corresponding  to  the  CFD-predicted  LEV  initiation 
(referred  to  as  uLESPrnaxv)  agree  with  the  critical  value  of  LESP  for  the  corresponding  2D  airfoil 
(referred  to  as  uLESPcrit”)7  And  how  does  the  spanwise  location  of  the  maximum  in  the  spanwise 
LESP  distribution  agree  with  the  CFD-predicted  spanwise  location  for  initiation  of  LEV  formation? 

(iii)  if  the  agreements  in  (ii)  are  good,  can  the  critical  value  of  LESP,  obtained  from  2D  studies,  be 
used  to  predict  the  time  instant,  angle  of  attack,  and  spanwise  location  for  LEV  initiation  using 
UVLM  (without  using  3D  CFD  studies),  and  how  closely  will  these  predicted  values  agree  with  the 
corresponding  CFD  predictions? 

The  following  section  presents  the  methodology  for  carrying  out  the  LEV  initiation  study  and  the 
case  studies  used  to  assess  the  effectiveness  of  the  low-order  prediction  using  the  UVLM  code.  The 
subsequent  section  presents  the  results  from  the  low-order  and  high-order  methods. 


8.1  Methodology  and  Case  Studies 

A  total  of  13  finite- wing  geometries  and  two  airfoil  sections  are  considered  in  this  effort.  The  two 
airfoil  sections  are  the  SD  7003  (see  the  reference107)  and  a  modified  SD  7003  with  a  sharpened 
leading  edge  having  a  50%  reduction  in  the  leading-edge  radius  compared  to  the  original  SD  7003 
airfoil.  The  two  airfoil  sections,  referred  to  as  the  SD7003  and  sharpened  SD7003  in  this  thesis,  are 
shown  in  figure  8.1. 

The  13  finite  wing  cases,  labelled  cases  1  13,  have  different  taper  ratios,  tip-twist  angles,  aspect 
ratios,  and  pivot  locations,  with  sections  formed  using  one  or  both  of  the  two  airfoils,  SD7003  and 
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(a)  Leading  edge. 
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(b)  Overview. 


Figure  8.1:  Comparison  of  SD  7003  (solid  line)  and  sharpened  SD  7003  (dashed  line) 


sharpened  SD7003.  Table  8.1  lists  the  details  of  the  two  airfoils  and  the  13  finite- wing  cases  used 
in  this  paper.  Figures  8.2  show  the  ten  wing  geometries  that  are  used  in  the  13  cases. 

Table  8.1:  Test  cases. 


Case 

Taper 

AR 

Pitch 

Pivot 

Tip 

Sweep 

Inboard 

Airfoil 

ratio 

rate 

point 

twist 

angle 

incidence 

K 

(x/cave) 

[deg] 

[deg] 

[deg] 

2D1 

- 

2D 

0.3 

0.25 

- 

- 

0 

SD7003 

2D2 

- 

2D 

0.3 

0.25 

- 

- 

0 

Sharpened  SD7003 

1 

1 

6 

0.3 

0.25 

0 

0 

0 

SD7003 

2 

1 

6 

0.3 

0.75 

0 

0 

0 

SD7003 

3 

1 

6 

0.2 

0.25 

0 

0 

0 

SD7003 

4 

1 

6 

0.4 

0.25 

0 

0 

0 

SD7003 

5 

0.5 

6 

0.3 

0.25 

0 

0 

0 

SD7003 

6 

1 

6 

0.3 

0.25 

10 

0 

0 

SD7003 

7 

1 

2 

0.3 

0.25 

0 

0 

0 

SD7003 

8 

1 

4 

0.3 

0.25 

0 

0 

0 

SD7003 

9 

1 

8 

0.3 

0.25 

0 

0 

0 

SD7003 

10 

1 

6 

0.3 

0.25 

0 

30 

0 

SD7003 

11 

1 

6 

0.3 

0.25 

0 

0 

0 

Sharpened  SD7003 

12 

1 

6 

0.3 

0.25 

0 

0 

4a 

SD7003 

13 

1 

6 

0.3 

0.25 

0 

0 

0 

SD7003" 

a)  Inboard  third  of  wing  has  a  4-degree  larger  incidence  compared  to  the  rest  of  the  wing. 

b)  Inboard  third  of  wing  has  the  sharpened  SD7003  airfoil,  with  the  SD7003  used  on  the  rest  of  the  wing. 


All  the  studies  in  this  work  have  been  performed  for  a  chord  Reynolds  number  of  20,000.  This 
value  was  chosen  because  our  previous  studies  have  shown  that  the  LESP  criterion  successfully 
predicts  initiation  of  LEV  formation  on  airfoils  (in  two-dimensional  flow)  at  Reynolds  numbers 
between  10,000  and  40,000.  In  this  range  of  Reynolds  numbers,  the  RANS  CFD  analysis  using 
the  Spalart-Allmaras  turbulence  model,  as  implemented  in  REACTMB-INS  flow  solver,  has  also 
been  shown  to  agree  well  with  experimental  results  for  LEV  initiation  and  formation  on  airfoils. 
The  only  case  with  a  non-constant  chord  is  case  5;  for  this  case,  the  leading-edge  has  zero  sweep, 
and  the  average  chord  (at  the  mid-semi-span  location)  has  been  used  as  the  length  scale  to  set  the 
Reynolds  number  and  the  non-dimensional  pitch  rate.  For  the  tapered  wing  (case  5),  the  pivot 


74 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


(a)  Rectangular  wing  (AR  2)  for  case  7. 


(c)  Rectangular  wing  (AR  6)  for  cases  1,  2,  3, and 
4. 


(b)  Rectangular  wing  (AR  4)  for  case  8. 


(d)  Rectangular  wing  (AR  8)  for  case  9. 


(e)  Tapered  wing  for  case  5.  (f)  Twisted  wing  for  case  6. 


(g)  Swept  wing  for  case  10. 


(h)  Rectangular  wing  with  sharpened  leading 
edge  for  case  11. 


(i)  Rectangular  wing  with  inboard  third  having  (j)  Rectangular  wing  with  inboard  third  having 
larger  incidence  for  case  12.  sharpened  leading  edge  for  case  13. 


Figure  8.2:  Geometries  of  the  ten  wings  used  in  the  13  cases. 


75 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


location  for  the  pitching  motion  is  at  the  25%  of  the  chord  at  the  mid-semi-span  location,  while 
for  the  swept  wing  (case  10),  the  pivot  is  at  the  quarter-chord  location  of  the  root  section.  Cases 
12  and  13  comprise  wings  that  have  abrupt  changes  in  geometry  demarcating  the  inboard  third  of 
the  wing  span  from  the  outboard  regions.  Case  12  has  a  4-degree  higher  incidence  on  the  inboard 
third  of  the  wing  compared  to  the  rest  of  the  wing,  and  case  13  has  the  sharpened  SD7003  airfoil 
over  the  inboard  third  of  the  airfoil  and  the  original  SD7003  section  on  the  outboard  portions.  The 
swept- wing  geometry  in  case  10  has  been  defined  using  the  airfoil  section  parallel  to  the  plane  of 
symmetry. 


8.1.1  Motion  parameters 

Although  the  LESP  criterion  has  been  verified  for  arbitrary  pitching,  plunging,  surging,  and  com¬ 
bination  motions,2  the  current  work,  owing  to  the  large  number  of  geometry  cases,  focusses  on  a 
pure  pitching  motion.  Thus,  for  all  wing  geometries  in  this  work,  a  0-45-degree  pitch-ramp  motion 
is  considered,  with  a  non-dimensional  pitch  rate  of  K  =  0.3  used  in  most  cases,  except  for  cases  3 
and  4,  which  have  K  =  0.2  and  0.4,  respectively.  Figure  8.3  shows  the  time  variations  of  the  angle 
of  attack  (same  as  pitch  angle  in  this  work)  for  the  three  pitch  rates.  The  equation  for  the  pitch 
ramp  motion  is  from  Granlund  et  al. 10 


Figure  8.3:  Pitch-up  motions  used  in  this  study. 


8.1.2  Determination  of  LEV  initiation  and  LESPmax  from  CFD  and  UVLM 

An  important  element  of  the  current  work  is  the  quantitative  determination  of  the  time  instant  of 
LEV  initiation  from  CFD  results.  In  past  research,  experimental  studies  (108  and,38  for  example) 
have  used  the  movement  of  the  minimum-pressure  location  to  track  movement  of  the  LEV,  and 
the  computational  study  of  Ref.48  brought  to  light  the  behavior  of  critical  points  in  the  velocity 
field  near  an  LEV.  Guided  by  these  results  in  the  literature,  in  our  earlier  work  on  LEV  initiation 
on  airfoils,2  we  showed  that  a  skin- friction  signature  near  the  leading  edge  from  CFD  results  could 
be  consistently  used  to  identify  LEV  initiation.  In  the  current  work,  we  adapt  this  skin-friction 
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signature  to  predict  LEV  initiation  on  finite-wing  flows.  In  the  remainder  of  this  subsection,  we 
discuss  this  skin-friction  signature,  first  for  airfoil  flows  and  next  for  finite-wing  flows. 

In  order  to  provide  an  overview  of  the  events  leading  up  to  LEV  formation,  Figures  8.4  shows 
a  series  of  representative  CFD  snapshots  of  streamlines  and  upper-surface  Cf  plots  for  the  SD 
7003  undergoing  pitch  up  motion  (Case  2D1).  At  the  start  of  the  pitching  motion,  the  airfoil  has 
attached  flow  over  most  of  the  upper  surface,  with  the  upper-surface  Cf  becoming  negative  past 
x/c  =  0.7  indicating  the  presence  of  a  small  region  of  reversed  flow  over  the  aft  30%  of  the  chord 
in  Figure  8.4a.  At  the  higher  angle  of  attack  of  18.2  degrees,  as  seen  from  Figure  8.4b,  the  trailing- 
edge  reversed  flow  region  extends  from  x/c  —  0.5.  Of  interest,  however,  is  the  tiny  region  near 
the  leading  edge  over  which  Cf  is  negative,  indicating  the  beginning  of  flow  reversal  at  the  leading 
edge.  At  a  higher  pitch  angle  of  23.4  degrees,  the  Cf  distribution  in  Figure  8.4c  shows  a  positive 
spike  reaching  up  to  Cf  =  0  within  the  negative-C/  region  near  the  leading  edge.  In  the  approach 
developed  in  our  earlier  work,2  this  first  occurrence  of  positive  Cf  within  the  negative-C/  region 
near  the  leading  edge  is  taken  as  the  time  instant  corresponding  to  initiation  of  LEV  formation. 
The  LEV  becomes  discernable  in  the  streamline  plot  at  the  higher  pitch  angle  of  26.0  degrees  in 
Figure  8.4d,  and  clearly  visible  at  even  higher  pitch  angles  (not  shown).  As  the  LEV  grows,  multiple 
vortices  near  the  primary  vortex  are  formed,  resulting  in  the  occurrence  of  several  positive  spikes 
with  the  negative-C/  region  near  the  leading  edge.  The  structures  observed  in  our  CFD  results  are 
similar  to  those  observed  by  Ghosh-Choudhuri  et  a/.48  In  the  current  work,  however,  the  focus  is 
on  the  initiation  of  LEV  formation  rather  than  on  the  flow  features  that  occur  subsequent  to  LEV 
initiation. 

In  extending  the  skin- friction  signature  for  LEV  initiation  to  finite- wing  flows,  we  examine  the  CFD 
plots  of  the  skin- friction  lines  on  the  upper-surface  at  successive  time  instants.  The  objective  is  to 
find  the  time  instant  corresponding  to  the  first  occurrence  of  a  region  of  positive  skin-friction  within 
the  negative  skin-friction  region  near  the  leading  edge.  To  illustrate  the  procedure,  Figures  8.5a 
and  8.5b  show  the  upper-surface  skin- friction  lines  on  the  wing  used  in  case  1  at  two  successive 
time  instants  from  CFD  output  corresponding  to  just  prior  to  LEV  initiation  and  just  after  LEV 
initiation,  respectively.  Figure  8.5a  indicates  that  there  are  roughly  four  flow  regions  at  t*  =  1.695: 
(i)  a  small  region  of  reversed  flow  near  the  leading  edge  with  negative  chordwise  Cf,  which  is  usually 
a  precursor  to  LEV  formation;  (ii)  a  thin  layer  of  reversed  flow  in  the  vicinity  of  the  trailing  edge, 
indicating  trailing-edge  flow  reversal;  (iii)  the  triangle-shaped  region  at  the  right  edge  resulting 
from  surface  flow  caused  by  the  tip  vortex;  and  (iv)  the  intermediate  flow  region  with  flow  having  a 
chordwise  component  that  is  from  leading  to  trailing  edge  corresponding  to  positive  chordwise  Cf. 

Figure  8.5b  shows  the  surface  streamlines  for  the  very  next  time  instant  (t*  =  1.710)  from  CFD 
output  for  case  1.  It  is  seen  that  there  is  a  new  region  near  the  leading  edge  of  the  root  area. 
This  small  region  is  the  first  occurrence  of  positive  chordwise  skin  friction  within  the  reversed-flow 
region  near  the  leading  edge.  By  analogy  to  the  skin-friction  signature  in  the  airfoil  case,  it  can 
be  said  that  the  t*  corresponding  to  LEV  initiation  for  the  finite- wing  case  1  is  between  1.695  and 
1.710.  It  is  also  seen  that  the  LEV  initiation  occurs  at  the  wing  root  for  this  wing,  i.e. ,  at  2 y/b  =  0, 
with  the  LEV  starting  to  form  over  a  spanwise  region  extending  approximately  from  2 y/b  =  —0.3 
to  0.3. 

Figure  8.6  shows  the  spanwise  distributions  of  LESP  for  case  1  for  t*  =  1.710  (solid  curve,  for 
CFD  frame  just  after  LEV  initiation)  and  for  t*  =  1.695  (dashed  curve,  for  CFD  frame  just  prior 
to  LEV  initiation).  These  spanwise  LESP  distributions  were  obtained  from  the  UVLM  analysis. 
The  spanwise-maximum  LESP  value,  corresponding  to  LEV  initiation,  for  case  1,  therefore,  occurs 
between  the  maximum  values  for  these  two  spanwise  LESP  distributions.  Thus,  the  solid  black 
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(a)  From  CFD  frame  just  prior  to  LEV  on-  (b)  From  CFD  frame  just  after  to  LEV  on¬ 
set,  t*  =  1.695,  a  =  23.88  deg.  set,  t*  =  1.710,  a  =  24.41  deg. 


Figure  8.5:  Upper-surface  skin-friction  lines  for  case  1.  Right  half  of  wing  shown.  In  each  snapshot, 
leading  edge  is  on  the  top,  wing  tip  is  on  the  right,  trailing  edge  is  on  the  bottom,  wing  root  is  on 
the  left. 


Nondimensional  spanwise  coordinate 


Figure  8.6:  Spanwise  variation  of  LESP  and  determination  of  LESPmax  for  case  1. 
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horizontal  line  in  figure  8.6  indicates  the  LESPmax  value,  with  the  black  error  bar  denoting  un¬ 
certainty,  corresponding  to  initiation  of  LEV  for  case  1.  Also  plotted  as  a  horizontal  line  is  the 
LESPcrit  from  2D  airfoil  CFD  result,  intended  to  be  compared  with  the  LESPmax  line. 

In  the  remainder  of  the  paper,  results  for  initiation  of  LEV  formation  on  any  arbitrary  wing  will  be 
presented  using  a  plot  of  upper-surface  skin-friction  lines  from  a  CFD  frame  corresponding  to  the 
instant  of  time  just  after  LEV  onset.  The  LESP  distribution  from  UVLM  will  also  be  presented  for 
the  t*  corresponding  to  this  frame.  The  LESPmax  for  the  wing  and  2D  LESPcru  for  the  airfoil  with 
error  estimates  will  be  denoted  by  horizontal  line  with  vertical  error  bars  in  the  LESP  vs.  spanwise 
coordinate  plots. 


8.2  Results  and  Discussion 

The  following  subsection  presents  the  values  of  LESPcra  for  the  SD  7003  (case  2D1)  and  the 
sharpened  SD  7003  (case  2D2)  airfoils  used  in  this  study.  The  next  subsection  presents  the  results 
for  the  baseline  finite- wing  geometry  (case  1),  in  which  the  finite- wing  results  are  compared  with 
the  2D  results  for  LEV  initiation.  The  subsequent  results  are  presented  as  seven  case  studies  to 
illustrate  the  effects  of  pivot  location  (case  study  A),  pitch  rate  (case  study  B),  wing  taper  ratio 
(case  study  C),  wing-tip  twist  (case  study  D),  aspect  ratio  (case  study  E),  sweep  angle  (case  study 
F),  leading-edge  shape  (case  study  G),  abrupt  change  in  section  incidence  (case  study  H),  and 
abrupt  change  in  section  leading-edge  radius  (case  study  I),  on  initiation  of  the  LEV  formation.  A 
summary  of  all  case  studies  is  then  discussed.  In  each  case  study,  case  1  is  used  as  the  baseline 
case  for  reference,  and  all  the  finite-wing  results  are  compared  with  the  corresponding  2D  results 
for  LEV  initiation. 

8.2.1  Two-dimensional  cases  2D1  and  2D2 

The  two  airfoils  used  in  the  study,  the  original  and  the  sharper-leading-edge  versions  of  the  SD 
7003,  were  studied  for  the  pitching  motions  listed  under  cases  2D1  and  2D2  in  Table  8.1.  For  each 
case,  results  from  two-dimensional  CFD  analysis  were  studied  to  determine  the  time  instant  and 
angle  of  attack  for  LEV  initiation  using  the  approach  described  in  Section  8.1.2.  The  time- variation 
of  LESP  for  each  case  was  determined  using  the  unsteady  thin  airfoil  theory  of  Ref.1  From  the 
results  for  case  2D1,  the  time  instant,  angle  of  attack  and  LESP  at  LEV  initiation  were  found  to 
be  1.680,  23.38  degrees,  and  0.269,  respectively.  Similarly,  the  results  for  case  2D2  yield  the  time 
instant,  angle  of  attack,  and  LESP  to  be  1.605,  20.80  degrees,  and  0.237  at  LEV  initiation.  Because 
case  2  uses  the  sharpened- leading-edge  airfoil,  the  LEV  initiation  occurs  at  an  earlier  time  in  the 
motion.  Thus,  the  LESPcra  values  for  the  SD7003  and  the  sharpened  SD7003  airfoils,  used  in  the 
remainder  of  this  report,  are  0.269  and  0.237. 

8.2.2  Baseline:  case  1 

The  results  for  case  1  are  presented  in  Figure  8.7.  Figure  8.7a,  shown  earlier  as  Figure  8.5b  and 
repeated  here  for  completeness,  shows  the  upper-surface  skin-friction  lines  from  CFD  analysis  at 
the  time  instant  just  after  LEV  initiation.  Using  the  skin-friction  signature  described  earlier  in 
Section.  8.1.2,  it  is  seen  that  the  LEV  starts  forming  at  the  wing  root.  Figures  8.7b,  8.7c,  and 
8.7d  show  the  leading-edge  flow  streamlines  (for  velocities  relative  to  the  body  frame)  at  the  tip, 
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(a)  Upper-surface  skin-friction  lines  for  case  1, 
t*  =  1.710,  a  =  24.41  deg. 


(b)  Streamlines  near  wing  tip,  leading  edge, 


(c)  Streamlines  at  2y/b  =  0.5,  leading  edge. 


(d)  Streamlines  at  wing  root,  leading  edge. 


Nondimensional  spanwise  coordinate 

(e)  LESP  distribution  at  LEV  initiation. 
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(g)  LESPmax  vs.  t*  from  UVLM. 


(h)  CL  vs.  t*  from  CFD  and  UVLM. 


Figure  8.7:  Result  of  baseline  (case  1) 
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mid-semi-span,  and  root  locations,  respectively.  These  streamline  plots  show  the  LEV  formation, 
in  the  early  stages,  visible  only  over  the  inboard  portion  of  the  wing.  The  fact  that  the  surface 
skin- friction  plot  in  Figure  8.7  does  not  show  evidence  of  LEV  formation  at  2 y/b  —  0.5  is  probably 
because  of  a  small  time  lag  in  the  signature  registering  on  the  surface. 

The  general  features  of  LEV  structure  at  the  initiation  for  the  baseline  case  agree  with  the  obser¬ 
vations  of  Freymuth,57  Yilmaz  and  Rockwell59  and  Spentzos  et  a/.109,110  in  that,  for  a  pitching 
rectangular  wing,  the  LEV  structure  is  parallel  to  the  leading  edge  over  the  inboard  portion  of  the 
wing,  with  its  maximum  strength  around  wing  root.  The  corner  portion  of  LEV  structure  tends 
to  attach  to  the  wing  surface,  whereas  the  inboard  portion  of  LEV  structure  is  lifted  up  and  gets 
convected  downstream,  as  reported  by  Granlund  et  al111  As  shown  in  figures  3.5,  the  LEV  subse¬ 
quently  forms  a  three-dimensionally  distorted  vortex  structure,  which  is  often  called  an  O-shaped 
structure,  bounded  at  the  leading  edge  of  the  wing  tips.  This  characteristic  LEV  structure  has 
been  observed  in  different  motion  kinematics  such  as  plunging  maneuver  by  Visbal  et  al.112  and 
the  flapping  cycle  maneuver  by  Gordnier  and  Demasi.113 

Figure  8.7e  shows  the  spanwise  distribution  of  LESP  on  the  right  side  of  the  wing  at  the  time 
instant  of  LEV  initiation.  The  occurrence  of  maximum  LESP  at  the  wing  root  agrees  well  with 
the  initiation  of  LEV  at  the  wing  root  in  the  CFD  prediction.  This  spanwise  LESP  distribution 
is  also  consistent  with  the  variation  of  suction  pressure  on  leading  edge  observed  by  Schreck  and 
Helin.58  Also  plotted  are  the  horizontal  line  at  the  maximum  value  and  an  error  bar  denoting 
the  change  in  maximum  value  between  two  consecutive  output  frames  from  CFD.  The  value  of 
LESPcrit  for  the  2D  airfoil  is  also  plotted  using  a  dashed  line  with  its  error  bar.  The  excellent 
agreement  between  the  3D  LESPmax  and  2D  LESPcrit  shows  that  the  initiation  of  LEV  formation 
on  the  finite  wing  could  have  been  predicted  using  a  low-order  analysis  and  the  LESPcru  value. 
This  type  of  low-order  prediction  is  referred  to  as  “UVLM  prediction”  in  the  remainder  of  this 
chapter.  Figure  8.7f  shows  the  motion  history  for  this  case  with  comparison  of  the  prediction  of 
LEV  initiation  from  UVLM  and  CFD.  The  agreement  is  seen  to  be  excellent. 

The  time  variation  of  LESPmax  from  UVLM  is  shown  in  Figure  8.7g,  with  the  dotted  vertical 
line  denoting  the  time  instant  of  LEV  initiation  from  CFD  prediction,  and  dotted  horizontal  line 
representing  the  LESPcrit  from  2D  CFD.  It  is  seen  that  LESPmax  increases  with  constant  slope 
during  pitch-ramp  motion.  Figure  8.7h  shows  the  lift-coefficient  variation  with  time.  At  the  start 
of  pitching,  both  CFD  and  UVLM  predict  a  spike  in  Cl-  This  is  due  to  apparent-mass  effects. 
During  the  linear  pitch  up,  although  Cl  from  CFD  and  UVLM  monotonically  increase  with  the 
angle  of  attack,  the  CFD  result  is  lower  than  UVLM  because  UVLM  is  based  on  inviscid  flow  while 
CFD  is  a  viscous-flow  solution.  However,  after  the  onset  of  LEV,  the  Cl  from  CFD  is  higher  than 
that  from  UVLM.  After  LEV  initiation,  CFD  takes  into  consideration  the  vortex  lift  due  to  LEV 
formation,  while  current  UVLM  does  not  model  the  LEV  effects.  Overall,  it  is  seen  that  prior  to 
LEV  initiation,  the  Cl  prediction  between  CFD  and  UVLM  are  in  good  agreement. 


8.2.3  Case  Study  A:  Effect  of  pivot  location 

In  case  study  A,  the  initiation  of  the  LEV  for  the  baseline  case  (case  1,  pivot  at  x/c  =  0.25)  is 
compared  with  that  for  the  same  wing  with  pivot  at  x/c  =  0.75  (case  2). 

Figure  8.8  shows  the  comparison  of  results  for  the  two  pivot  locations.  While  the  spanwise  location 
for  initiation  of  the  LEV  for  case  2  (as  deduced  from  the  skin-friction  plot  in  Figure  8.8a)  is  similar 
to  that  for  the  baseline  case  1,  it  is  seen  that  case  2  has  a  larger  region  of  trailing-edge  reversed 
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(a)  Upper-surface  skin-friction  lines  for  case  2,  t*  =  2.205,  a  =  41.14  deg. 


Nondimensional  spanwise  coordinate 


(b)  LESP  distribution  at  LEV  initiation. 
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(d)  LESPmax  vs.  t*  from  UVLM. 


Figure  8.8:  Case  study  A:  Effect  of  pivot  location.  Comparison  of  cases  1  and  2  at  CFD  frames 
just  after  LEV  onset. 
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flow  that  extends  to  almost  the  reversed  flow  at  the  leading  edge.  The  major  difference  from  the 
baseline,  however,  is  that  LEV  initiation  is  largely  delayed  due  to  motion-induced  “downwash” 
at  the  leading  edge.  This  trend  of  delayed  LEV  formation  with  aft  movement  of  the  pivot  point 
is  not  new.  The  trend  and  the  reasoning  have  been  reported  by  several  researchers  including 
Ham  and  Garelick,45  Visbal  and  Shang,46  Visbal  and  Gordnier,114  01, 49  and  Granlund  et  a/.10,111 
Furthermore,  according  to  Ericsson,115  rearward  pivot  location  increases  the  pitching  speed  at  the 
leading  edge  and  the  moving  wall  effect  gives  more  beneficial  effect  in  boundary  layer  around  the 
leading  edge  and  thus  delays  LEV  formation. 

In  figure  8.8b,  the  spanwise  variation  of  LESP  is  similar  to  that  for  the  baseline,  but  there  is  slight 
discrepancy  between  the  LESPmax  for  case  2  and  the  2D  LESPcrit .  This  discrepancy  manifests  as 
an  approximately  1.2-degree  discrepancy  between  the  CFD  and  UVLM  predictions  for  the  AoA  for 
LEV  initiation  (Figure  8.8c).  The  motion-induced  downwash  due  to  the  aft  pivot  location  in  case  2 
also  reduce  the  LESPmax  (Figure  8.8d)  and  the  lift  coefficient  (Figure  8.8e)  in  comparison  to  those 
for  case  1.  As  seen  in  case  1,  the  UVLM-predicted  Cl  for  case  2  is  higher  than  the  CFD-predicted 
value  until  LEV  initiation,  beyond  which  the  UVLM  predictions  become  progressively  invalid. 

8.2.4  Case  Study  B:  Effect  of  pitch  rate 

In  case  study  B,  the  initiation  of  the  LEV  for  the  baseline  case  (case  1,  K  =  0.3)  is  compared  with 
that  for  the  same  wing  with  K  =  0.2  (case  3)  and  K  =  0.4  (case  4). 

The  results  for  the  case  study  for  three  pitch  rates  is  shown  in  figure  8.9.  The  upper-surface  skin- 
friction  lines,  plotted  in  Figures  8.9a  and  8.9b  for  K  =  0.2  and  0.4,  respectively,  show  that  the  LEV 
initiation  starts  at  the  wing  root  like  for  case  1.  It  is  also  seen  that,  with  increasing  pitch  rate,  the 
chordwise  extent  of  the  tr ailing-edge  reversed- flow  region  progressively  reduces.  There  is  a  small 
increase  in  AoA  for  LEV  initiation  with  increase  in  pitch  rate.  Increase  in  pitch  rate  increases  the 
motion-induced  downwash  at  the  leading  edge,  causing  the  small  delay  in  LEV  initiation.  This 
trend  between  the  pitch  rate  and  the  time  instant  of  LEV  onset  is  consistent  with  the  features  and 
explanations  in  the  literature:  Ham  and  Garelick,45  Johnson  and  Ham,116  McCroskey  et  a/.,117  119 
Walker  et  a/.,120  Carr,121  Visbal  and  Shang,46  Acharya  and  Metwally,122  Schreck  and  Helin,58 
Choudhuri  and  Knight,47  Lian,123  and  Jantzen  et  al.124: 

Figure  8.9c  shows  the  spanwise  LESP  distributions  at  LEV  initiation  for  cases  3  and  4.  The 
maximum  values  for  the  LESP  distributions  are  at  2 y/b  =0,  which  agrees  with  the  CFD-predicted 
behavior  of  LEV  initiation  starting  at  the  wing  root.  The  figure  also  shows  that  the  LESPmax 
value  for  the  slow  pitch  rate  (case  3,  dashed  line)  is  lower  than  LESPcra  and  LESPmax  for  the 
fast  pitch  rate  (case  4,  dotted  and  dashed  line)  is  higher  than  LESPcrit.  As  a  result,  as  seen  in 
Figure  8.9d,  the  UVLM  predictions  for  the  AoA  for  LEV  initiation  for  cases  3  and  4  are  respectively 
higher  and  lower  than  the  CFD-predicted  values.  Figure  8.9e  shows  the  variations  of  LESPmax  for 
cases  1,  3,  and  4.  Figure  8.9f  shows  the  effect  of  pitch  rate  in  Cl-  Generally  speaking,  the  faster 
pitch  rate  generates  higher  Cl  and  this  trend  is  consistent  with  the  observations  of  McCroskey  and 
Pucci119  and  Granlund  et  a/.10,111 

8.2.5  Case  Study  C:  Effect  of  taper  ratio 

In  case  study  C,  the  initiation  of  the  LEV  for  the  tapered  wing  of  case  5  (taper  ratio  of  0.5,  with 
unswept  leading  edge)  is  compared  with  that  for  the  baseline  case. 
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(a)  Upper-surface  skin- friction  lines  for  case  3,  t*  =  1.980,  a  =  22.46  deg. 


(b)  Upper-surface  skin-friction  lines  for  case  4,  t*  =  1.562,  a  =  25.76  deg. 


Nondimensional  spanwise  coordinate 


(c)  LESP  distribution  at  LEV  initiation. 
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(e)  LESPmax  vs.  t*  from  UVLM. 


(f)  CL  vs.  t*  from  CFD  and  UVLM. 


Figure  8.9:  Case  study  B:  Effect  of  pitch  rate.  Comparison  of  cases  1,  3,  and  4  at  CFD  frames  just 
after  LEV  onset. 
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(a)  Upper-surface  skin-friction  lines  for  case  5, 
t*  =  1.725,  a  =  24.92  deg. 


Nondimensional  spanwise  coordinate 

(e)  LESP  distribution  at  LEV  initiation. 


(b)  Streamlines  near  wing  tip,  leading  edge. 


(c)  Streamlines  at  2y/b  =  0.3,  leading  edge. 


(d)  Streamlines  at  wing  root,  leading  edge. 


(f)  Motion  histories  and  LEV  initiation. 
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(g)  LESPmax  vs.  t*  from  UVLM. 


(h)  CL  vs.  t*  from  CFD  and  UVLM. 


Figure  8.10:  Case  study  C:  Effect  of  taper  ratio.  Comparison  of  cases  1  and  5  at  CFD  frames  just 
after  LEV  onset. 
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Figure  8.10a  shows  the  skin- friction  lines  on  upper  surface  for  case  5  at  LEV  initiation.  It  is 
seen  that  the  LEV  initiation  starts  at  around  2 y/b  —  0.3,  which  is  qualitatively  confirmed  by  the 
streamline  plots  around  the  leading  edge  at  the  three  spanwise  locations  in  Figures  8.10b  8.  lOd  and 
is  similar  to  the  observation  of  Spentzos  et  al110  that  LEV  initiation  on  a  tapered  wing  appeared 
to  be  closer  to  the  wing  tip.  Figure  8.10e  compares  the  spanwise  LESP  distribution  for  case  5 
with  that  for  case  1.  In  contrast  to  the  results  for  case  1,  the  spanwise  LESP  distribution  for  the 
tapered  wing  has  a  maximum  near  2 y/b  =  0.3,  which  agrees  well  with  the  CFD-predicted  spanwise 
location  for  LEV  initiation.  The  LESPmax  value  also  agrees  well  with  the  2D  LESPcru  value. 
Because  of  this  agreement,  the  UVLM-predicted  AoA  for  LEV  initiation  agrees  excellently  with 
the  CFD  prediction  (Figure  8.10f).  Figures  8.10g  and  8.10h  show  that  the  time  variations  of  the 
LESPmax  and  Cl  for  case  5  are  largely  similar  to  those  for  case  1. 

8.2.6  Case  Study  D:  Effect  of  tip  twist 

In  case  study  D,  the  initiation  of  the  LEV  for  the  twisted  wing  with  10-degree  tip  twist  (case  6)  is 
compared  with  the  results  for  case  1.  The  twisted  wing  in  case  6  has  a  linearly  increasing  section 
incidence  angle  from  root  to  tip. 

Figure  8.11a  shows  the  upper-surface  skin- friction  lines  from  CFD  at  the  time  instant  of  LEV 
initiation,  showing  that  LEV  starts  forming  for  this  wing  from  2 y/b  of  approximately  0.6.  The 
chordwise  streamline  plots  for  the  three  sections,  shown  in  Figures  8.11b,  8.11c,  and  8. lid,  confirm 
this  observation.  As  seen  in  Figure  8.11g,  the  UVLM-predicted  spanwise  LESP  distribution  of 
case  6  (dashed  line)  has  a  clear  maximum  around  the  spanwise  location  of  2 y/b  of  0.6,  which  agrees 
excellently  with  the  CFD  prediction.  The  LESPmax  value  from  UVLM  agrees  excellently  with 
the  2D  LESP crn  value.  As  a  result  of  this  agreement,  the  UVLM-predicted  value  for  the  AoA  for 
LEV  initiation  also  agrees  excellently  with  the  CFD  prediction,  as  seen  from  Figure  8. Ilf.  From 
Figure  8. Ilf,  it  is  also  seen  that  the  AoA  for  LEV  initiation  for  case  6  is  approximately  five  degrees 
less  than  that  for  case  1.  This  early  initiation  of  the  LEV  formation  is  the  result  of  higher  incidence 
of  the  outboard  sections  due  to  tip  twist.  The  behavior  of  the  LESPmax  and  Cl  time  variations, 
shown  in  Figures.  8.11g  and  8.11h,  are  largely  similar  to  those  for  case  1,  except  for  an  offset 
resulting  from  the  wing  twist. 

8.2.7  Case  Study  E:  Effect  of  aspect  ratio 

In  case  study  E,  the  initiation  of  the  LEV  for  the  rectangular  wing  cases  7,  8,  and  9  ( AR  =  2,  4,  8 
respectively)  are  compared  with  that  for  the  baseline  case  (case  1,  rectangular  wing  of  AR  =  6). 

The  results  of  the  investigation  of  AR  effects  on  LEV  initiation  are  presented  in  Figure  8.12. 
Figures  8.12a,  8.12b,  and  8.12c  show  the  skin  friction  lines  for  the  rectangular  wings  with  AR  = 
2,  4,  and  8,  respectively.  These  three  figures  indicate  that  the  onset  location  of  LEV  initiation  is 
at  the  wing  root  for  all  these  wings.  Also  seen  is  that,  although  the  region  near  the  wing  tip  for 
each  wing  where  the  skin-friction  lines  are  influenced  by  the  tip  vortex  are  roughly  of  the  same  size 
for  these  three  wings,  the  fraction  of  the  wing  occupied  by  this  region  is  clearly  the  largest  for  the 
lower  aspect  ratio  wing.  Thus  it  can  be  expected  that  the  lower  aspect  ratio  wings  will  be  more 
influenced  by  the  tip  vortex  effects.  Figure  8.12d  shows  the  spanwise  distributions  of  LESP  for 
the  four  cases:  AR  =  2,  4,  6,  and  8.  Each  of  the  four  cases  has  the  maximum  in  LESP  variation 
at  the  root,  which  agrees  with  the  CFD-predicted  root  location  for  the  onset  of  LEV  formation  for 
each  wing.  While  the  AR  4,  6,  and  8  wings  have  LESPmax  values  agreeing  excellently  with  the  2D 
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(a)  Upper-surface  skin-friction  lines  for  case  6, 
t*  =  1.575,  a  =  19.75  deg. 


(b)  Streamlines  near  wing  tip,  leading  edge. 


(c)  Streamlines  at  2y/b  =  0.6,  leading  edge. 


(d)  Streamlines  at  wing  root,  leading  edge. 
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(e)  LESP  distribution  at  LEV  initiation. 


(f)  Motion  histories  and  LEV  initiation. 
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(g)  LESPmax  vs.  t*  from  UVLM. 


(h)  CL  vs.  t*  from  CFD  and  UVLM. 


Figure  8.11:  Case  study  D:  Effect  of  tip  twist.  Comparison  of  cases  1  and  6  at  CFD  frames  just 
after  LEV  onset. 
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(a)  Upper-surface  skin- friction  lines  for  case  7,  t* 


1.830,  a 


(b)  Upper-surface  skin- friction  lines  for  case  8,  t*  =  1.740,  a 
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(c)  Upper-surface  skin- friction  lines 


for  case  9,  t*  =  1.695,  a  =  23.89  deg. 
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(d)  LESP  distribution  at  LEV  initiation. 


(e)  Motion  histories  and  LEV  initiation. 
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(f)  LESPmax  vs.  t*  from  UVLM. 


(g)  CL  vs.  t*  from  CFD  and  UVLM. 


Figure  8.12:  Case  study  E:  Effect  of  aspect  ratio.  Comparison  of  cases  1,  7,  8,  and  9  at  CFD  frames 
just  after  LEV  onset. 
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(c)  LESPmax  VS.  t*  from  UVLM. 


(d)  CL  vs.  F  from  CFD  and  UVLM. 


Figure  8.13:  Comparison  of  separated  tip  flow  model  and  attached  tip  flow  model  for  case  1  and  7. 


LESPcrit,  the  LESPmax  for  the  AR  2  wing  is  noticeably  less  than  the  2D  LESPcra.  In  concert 
with  this  observation,  Figure  8.12e  shows  that  the  AoA  for  LEV  initiation  for  the  AR  4,  6,  and 
8  wings  from  both  CFD  and  UVLM  predictions  are  very  close  to  each  other.  On  the  other  hand, 
the  AR  2  wing  has  CFD  prediction  for  LEV  initiation  occurring  at  a  higher  AoA  compared  to  the 
UVLM  prediction.  This  discrepancy  can  be  attributed  to  the  increased  effect  of  the  “lifted-up”  tip 
vortex  structure  on  the  downwash  at  the  leading  edge  for  the  AR-2  case  due  to  the  proximity  of 
the  tip  vortex  to  the  wing  root.  Because  the  UVLM  does  not  take  the  lifting- up  of  the  tip  vortex 
into  consideration,  and  instead  assumes  an  attached  tip  vortex,  the  UVLM-predicted  AoA  for  LEV 
initiation  on  the  AR-2  wing  differs  from  the  CFD  prediction  by  almost  2.5  degrees.  Figures  8.12f 
and  8.12g  show  the  time  variations  of  LESPmax  and  Cl  for  the  four  cases.  Except  for  the  AR-2 
wing,  the  results  for  the  other  cases  mostly  follow  expected  trends. 

To  investigate  the  influence  of  the  tip- vortex  structure  on  the  UVLM  predictions,  the  UVLM  code 
was  used  to  study  the  AR-2  (case  7)  and  AR-6  (case  1)  wings  using  the  “attached  tip-flow  model, 
or  ATFM”  option  (Figure  7.4a)  and  the  “separated  tip-flow  model,  or  STFM”  option  (Figure  7.4b). 
Figure  8.13  presents  the  results  from  this  study.  It  is  seen  that  the  results  for  AR-6  wing  are  mostly 
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independent  of  the  tip-flow  model  used.  On  the  other  hand,  for  the  AR-2  wing,  the  results  of  the 
UVLM  have  better  agreement  with  the  CFD  prediction  when  the  separated  tip-flow  model  is  used 
than  when  the  attached  tip-flow  model  is  used.  This  trend  confirms  the  hypothesis  that,  for  very 
low  aspect  ratio  wings,  it  is  important  to  capture  the  effects  of  the  “lifted-off”  tip  vortex  structure 
for  improved  predictions  from  the  UVLM. 

8.2.8  Case  Study  F:  Effect  of  sweep  angle 

In  case  study  F,  the  initiation  of  the  LEV  for  the  30-degree  swept-wing  case  10  is  compared  with 
the  baseline  case  (case  1,  0-degree  sweep).  It  is  to  be  noted  that  the  pivot  point  for  the  swept  wing 
is  at  the  quarter-chord  location  of  the  root  section. 

The  results  for  the  effects  of  sweep  on  LEV  initiation  are  shown  in  figure  8.14.  Figure  8.14a  shows 
the  CFD-predicted  skin-friction  lines  on  the  upper-surface  at  LEV  initiation.  In  contrast  to  the 
results  for  all  the  other  wings,  this  swept  wing  case  exhibits  LEV  initiation  that  starts  close  to  the 
wing  tip,  and  at  a  very  small  AoA  of  2.8  degrees.  The  streamline  plots  of  the  leading-edge  regions 
in  Figures  8.14b,  8.14c,  and  8.14d  show  an  LEV  structure  for  the  sections  near  the  wing  tip  and 
none  for  the  root  section,  confirming  this  observation.  The  likely  reason  for  this  behavior  is  that  a 
pitch-up  motion  about  the  root  quarter-chord  location  causes  a  significant  motion-induced  upwash 
towards  the  wing  tips,  causing  the  leading  edge  near  the  wing  tip  to  become  critical  even  at  an  early 
stage  of  the  motion.  Figure  8.14e  compares  the  spanwise  distributions  of  LESP  at  LEV  initiation 
for  the  swept  wing  with  the  unswept  wing.  It  is  seen  that  the  LESP  distribution  for  the  swept 
wing  is  distinctly  different  from  that  for  the  unswept  wing;  it  has  a  clear  maximum  near  2 y/b  =  0.9, 
which  shows  that,  like  CFD,  the  UVLM  also  predicts  LEV  initiation  very  close  to  the  wing  tip. 
There  is  a  significant  discrepancy  between  the  LESPmax  for  the  swept  wing  and  the  2D  LESPcrit. 
One  possible  reason  for  this  discrepancy  is  that,  like  with  the  AR-2  wing,  the  LEV  occurs  very 
close  to  the  tip  vortex.  Because  the  UVLM  models  the  tip  vortex  using  an  “attached-tip-flow” 
model  (Figure  7.4a),  there  is  error  in  the  predicted  value  of  LESP  near  the  wing  tips.  Another 
reason  could  be  that  spanwise  pressure  gradient  or  spanwise  flow,  which  have  both  been  discussed 
by  other  researchers  in  the  context  of  LEV  initiation,33,60,125  could  be  causing  the  discrepancy. 

Figure  8.14f  compares  the  CFD-  and  UVLM-predicted  AoA  for  LEV  initiation  on  the  swept  wing 
with  those  for  the  unswept  wing  (case  1).  In  spite  of  the  discrepancy  between  the  LESPmax  and 
2D  LESPcrit  for  case  10 ,  it  is  seen  that  the  UVLM-predicted  AoA  is  a  small  value  of  1.5  degrees, 
which  is  within  1.3  degrees  of  the  CFD-predicted  AoA.  Figure  8.14g  shows  the  time- variation  of 
LESPmax  for  the  swept  and  unswept- wing  cases.  In  contrast  to  the  other  cases,  the  swept  wing 
is  seen  to  have  a  rapid  increase  in  the  LESPmax  even  during  the  very  early  stages  of  the  pitch- up 
motion.  Perhaps  for  this  reason,  the  fairly  significant  discrepancy  of  0.08  between  LESPmax  and 
2D  LESPcrit  translates  to  only  a  small  AoA  discrepancy  of  1.3  degrees  between  UVLM  and  CFD. 
Figure  8.14h  shows  the  time  variation  of  Cl  from  UVLM  and  CFD  for  the  two  cases.  Because 
the  LEV  initiation  occurs  very  early  in  the  motion  and  UVLM  does  not  model  the  LEV  formation 
and  its  effects,  there  is  considerable  lack  of  agreement  between  UVLM  and  CFD  for  much  of  the 
pitch-up  portion  of  the  motion. 

8.2.9  Case  Study  G:  Effect  of  leading-edge  curvature 

In  case  study  G,  the  initiation  of  the  LEV  for  the  wing  with  sharpened  leading  edge  (case  11, 
sharpened  SD7003  airfoil)  is  compared  with  the  baseline  case  (case  1,  SD7003). 
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(a)  Upper-surface  skin-friction  lines  for  case  10, 
t*  =  1.080,  a  =  2.79  deg. 


(b)  Streamlines  near  wing  tip,  leading  edge. 


(c)  Streamlines  at  2y/b  =  0.9,  leading  edge. 


(d)  Streamlines  at  wing  root,  leading  edge. 
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(e)  LESP  distribution  at  LEV  initiation.  (f)  Motion  histories  and  LEV  initiation. 
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(g)  LESPmax  vs.  t*  from  UVLM. 


(h)  CL  vs.  t*  from  CFD  and  UVLM. 


Figure  8.14:  Case  study  F:  Effect  of  sweep  angle.  Comparison  of  cases  1  and  10  at  CFD  frames 
just  after  LEV  onset. 
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(a)  Upper-surface  skin-friction  lines  for  case  11,  t*  =  1.620,  a  =  21.31  deg. 


Nondimensional  spanwise  coordinate 


(b)  LESP  distribution  at  LEV  initiation. 
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(d)  LESPmax  vs.  t*  from  UVLM. 


(e)  CL  vs.  £*  from  CFD  and  UVLM. 


Figure  8.15:  Case  study  G:  Effect  of  leading  edge  curvature.  Comparison  of  cases  1  and  11  at  CFD 
frames  just  after  LEV  onset. 
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Figure  8.15  presents  the  result  for  this  case  study.  The  results  largely  agree  with  the  expected 
behavior  that  airfoils  with  sharper  leading  edges  tend  to  have  earlier  leading-edge  separation  or  ini¬ 
tiation  of  vortex  shedding.39,44  Figure  8.15a  shows  the  upper-surface  skin- friction  lines  for  case  11. 
It  is  seen  that  LEV  initiation  starts  at  the  root,  like  with  case  1,  but  at  a  smaller  AoA  due  to  the 
sharper  leading  edge.  Compared  to  case  1,  the  chordwise  extend  of  trailing-edge  reverse  flow  region 
is  smaller  because  of  the  smaller  AoA.  The  spanwise  variation  of  LESP  for  the  two  cases  are  shown 
in  Figure  8.15b.  The  LESPmax  for  case  11  compared  excellently  with  the  2D  LESPcru  value  for  the 
sharpened  SD7003.  This  excellent  agreement  translates  to  excellent  agreement  between  the  UVLM 
and  CFD  predictions  for  the  AoA  for  LEV  initiation  for  case  11  (Figure  8.15c).  Figure  8.15d  shows 
that  the  time  variations  of  LESPmax  for  cases  1  and  11  are  identical.  Likewise  the  time  variations 
of  Cl  for  the  two  cases  from  UVLM  predictions  are  also  identical.  This  behavior  is  because  the  two 
airfoils  have  the  same  camberline,  and  the  UVLM  results  do  not  depend  on  leading-edge  radius. 
The  CFD  predictions  for  time  variation  of  Cl  for  case  11  is  seen  to  be  very  similar  to  that  for 
case  1. 


8.2.10  Case  Study  H:  Effect  of  abrupt  change  in  incidence 

In  case  study  H,  the  initiation  of  the  LEV  for  case  12,  for  which  the  wing  geometry  has  a  four-degree 
higher  incidence  over  the  inboard  third  of  the  span,  is  compared  with  that  for  case  1. 

Figure  8.16  shows  the  results  for  this  case  study.  The  upper-surface  skin- friction  lines  for  case  12  in 
Figure  8.16a  shows  the  abrupt  change  in  sectional  skin-friction  lines  at  2 y/b  =0.33,  which  is  clearly 
related  to  the  incidence  change.  It  is  seen  that  the  LEV  initiation  starts  at  the  wing  root,  which  is 
also  to  be  expected  owing  to  the  increased  incidence  in  that  region.  Figure  8.16b  shows  the  spanwise 
distributions  of  LESP  for  the  two  cases.  It  is  seen  that  case  12  has  a  distinctly  higher  LESP  over 
the  inboard  portion  of  the  wing.  The  LESPmax  for  case  12  is  seen  to  agree  closely  with  the  2D 
LESPcrit.  From  Figure  8.16c,  it  is  seen  that  the  UVLM  and  CFD  predictions  for  AoA  for  LEV 
initiation  on  case  12  agree  well  with  each  other,  and  that  case  12  has  LEV  initiation  approximately 
4  degrees  earlier  than  case  1.  Figure  8.16d  shows  that  the  time  variation  of  LESPmax  for  case  12  is 
quite  similar  to  that  of  case  1,  except  for  an  offset  resulting  from  the  increased  incidence.  Similarly, 
Figure  8.16e  shows  that  the  time  variation  of  Cl  for  case  12  is  similar  to  that  for  case  1,  except  for 
an  offset  due  to  the  incidence  change. 

8.2.11  Case  Study  I:  Partially  sharpened  wing 

In  case  study  I,  the  initiation  of  the  LEV  for  case  13,  with  the  wing  having  the  sharpened  leading- 
edge  over  the  inboard  third  of  the  wing  span,  is  compared  with  the  baseline  case  1. 

Figure  8.17  presents  the  results  for  this  case  study.  The  upper-surface  skin-friction  lines  for  case  13, 
in  Figure  8.17a,  shows  that  the  LEV  initiation  starts  at  the  wing  root  at  AoA  of  21.3  degrees. 
Figure  8.17b  compares  the  spanwise  variations  of  LESP  and  LESPcra  for  case  13  with  those  for 
the  baseline.  Because  of  the  abrupt  change  in  airfoil  between  the  inboard  and  the  outboard  portions 
of  the  wing,  the  2D  LESPcrit  value  also  changes  from  the  value  of  LESPcrn  =  0.269  for  the  SD 
7003  airfoil  in  the  outboard  portions  to  the  value  of  LESPcru  =  0.237  for  the  sharpened  SD  7003 
airfoil  in  the  inboard  portion.  It  is  seen  that  the  LESPmax  agrees  well  with  the  inboard  value  of 
the  2D  LESPcrit ,  resulting  in  an  excellent  agreement  between  the  UVLM  and  CFD  predictions 
for  the  AoA  for  LEV  initiation  (Figure  8.17c).  In  contrast  with  case  study  H  in  which  the  abrupt 
change  in  section  incidence  caused  a  sudden  change  in  the  spanwise  variation  of  LESP ,  the  abrupt 
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(a)  Upper-surface  skin-friction  lines  for  case  12,  t*  =  1.605,  a  =  20.80  deg. 


Nondimensional  spanwise  coordinate 


Nondimensional  time 


(b)  LESP  distribution  at  LEV  initiation.  (c)  Motion  histories  and  LEV  initiation. 
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(d)  LESPmax  vs.  t*  from  UVLM. 


(e)  CL  vs.  £*  from  CFD  and  UVLM. 


Figure  8.16:  Case  study  H:  Partially  deflected  wing.  Comparison  of  cases  1  and  12  at  CFD  frames 
just  after  LEV  onset. 
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Angle  of  attack  (AoA).  deg 


(a)  Upper-surface  skin-friction  lines  for  case  13,  t* 


1.620,  a  =  21.31  deg. 


Nondimensional  spanwise  coordinate 


(b)  LESP  distribution  at  LEV  initiation. 
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(d)  LESPmax  vs.  t*  from  UVLM. 


(e)  CL  vs.  t*  from  CFD  and  UVLM. 


Figure  8.17:  Case  study  I:  Partially  sharpened  wing.  Comparison  of  cases  1  and  13  at  CFD  frames 
just  after  LEV  onset. 
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(a)  Angle  of  attack  at  LEV  initiation. 


Figure  8.18:  Comparison  of  low-order  (UVLM) 
the  finite- wing  cases. 


High-order  prediction  of  spanwise  coordinate 

(b)  Spanwise  location  of  LEV  initiation, 
predictions  with  high-order  (CFD)  results  for  all 


change  in  leading-edge  radius  in  this  case  study  causes  an  abrupt  change  in  LESPcra.  Thus,  while 
both  cases  12  and  13  have  a  roughly  3-4-degree  earlier  LEV  initiation  compared  to  the  baseline 
case,  the  earlier  LEV  initiation  is  achieved  in  different  ways.  As  seen  from  Figure  8.17d,  the  time 
variations  in  LESPmax  for  both  cases  are  identical,  because  the  UVLM  prediction  is  not  altered 
by  changes  in  leading-edge  radius  as  long  as  section  camber  is  not  altered.  Figure  8.17e  shows  that 
the  time  variations  in  Cl  for  the  two  cases  are  also  nearly  identical. 

8.2.12  Summary 

Figure  8.18  presents  a  comparison  of  low-order  (UVLM)  predictions  for  all  the  13  finite- wing  cases 
with  the  respective  high-order  (CFD)  predictions.  In  Figure  8.18a,  the  AoA  values  for  LEV  ini¬ 
tiation  from  UVLM  are  plotted  as  symbols  for  the  13  cases  against  the  CFD  predictions.  The 
45-degree  line  is  co-plotted  to  denote  perfect  correlation,  with  the  thickness  of  the  line  chosen  to 
denote  the  measuring  error  in  the  CFD  observation.  Each  symbol  has  an  error  bar  to  denote  the 
resolution  error  in  the  UVLM  prediction.  It  is  seen  that  the  AoA  for  LEV  initiation  has  a  spread  of 
almost  40  degrees  varying  from  2.8  degrees  for  case  10  to  41.4  degrees  for  case  2.  This  shows  that 
AoA  itself  is  not  a  good  measure  for  LEV  initiation,  as  it  varies  significant  from  case  to  case.  In 
contrast,  when  the  low-order  prediction  is  made  using  the  2D  LESPcrit  value,  the  UVLM-to-CFD 
correlation  for  AoA  for  LEV  initiation  is  seen  to  be  very  good  for  all  the  cases  as  the  symbols  are 
all  close  to  the  45-degree  line.  The  error  in  UVLM-to-CFD  correlation  for  AoA  is  a  maximum  of 
2.6  degrees  for  case  7,  which  is  the  AR-2  wing,  with  most  of  the  other  cases  having  much  smaller 
error.  Figure  8.18b  shows  a  plot  of  the  low-order  (UVLM)  prediction  of  the  spanwise  location  of 
LEV  initiation  against  the  corresponding  CFD  prediction,  using  symbols  to  denote  the  13  cases. 
In  this  figure  too,  the  45-degree  line  is  plotted  to  show  perfect  correlation  with  the  thickness  of 
this  line  chosen  to  represent  the  resolution  error  of  the  UVLM.  The  spanwise  width  over  which 
LEV  initiation  was  observed  in  the  CFD,  averaged  over  the  13  cases,  is  shown  as  an  error  bar  in 
the  figure.  It  is  seen  that,  although  the  spread  for  the  spanwise  locations  ranges  from  the  wing 
root  to  almost  the  wing  tip,  the  UVLM  predictions  for  all  cases  agree  remarkably  well  with  the 
corresponding  CFD  observations.  These  excellent  correlations  in  (i)  AoA  for  LEV  initiation,  with  a 
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highest  error  of  2.6  degrees  for  cases  that  are  spread  over  a  range  of  41.4  degrees,  and  (ii)  spanwise 
location  with  a  highest  error  of  approximately  2 y/b  =  0.1  for  cases  that  are  spread  of  a  range  of 
2 y/b  =  0.9  demonstrate  the  power  of  the  LESP  concept  in  predicting  the  LEV  initiation  on  a  range 
of  finite- wing  geometries  using  just  the  2D  LESPcrit  values  for  the  sections  used  in  the  wings  along 
with  a  low-order  method  like  the  UYLM. 


8.3  Interim  Conclusions 

In  this  chapter  we  presented  a  study  to  extend  the  concept  of  LESP,  which  is  a  measure  of  leading 
edge  suction  in  2D  airfoils,  as  a  criterion  for  LEV  initiation  on  finite  wings.  In  order  to  assess 
the  effectiveness  of  this  idea,  the  maximum  values  of  the  spanwise  LESP  variations  for  the  wings 
at  LEV  initiation  were  studied  for  13  wings  with  variations  in  planform  geometry,  pivot  location, 
pitch  rate,  twist  angle,  aspect  ratio,  and  cross-section  shape.  It  was  seen  that,  for  any  given 
airfoil  shape  and  Reynolds  number,  this  maximum  LESP  value  for  all  the  wings  and  motions 
was  acceptably  close  to  the  critical  LESP  (obtained  from  2D  CFD).  It  was  shown  that  the  LESP 
concept  was  successful  in  predicting  the  angle  of  attack  for  the  onset  of  LEV  initiation  within  ±1.6 
degrees  accuracy  against  a  39-degree  spread  in  the  angle  of  attack  for  the  LEV  onset  for  the  cases 
examined.  Additionally,  the  LESP  concept  is  also  able  to  predict  the  spanwise  location  for  LEV 
initiation  with  good  accuracy.  Thus,  the  results  of  this  portion  of  the  research  effort  provided 
the  confidence  to  extend  the  LESP  concept  to  prediction  of  LEV  formation  (LEV  growth  after 
initiation),  which  is  discussed  in  the  next  chapter. 
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Chapter  9 


Low-Order  Prediction  of  LEV 
Formation  on  Finite  Wings 


This  chapter  presents  the  results  from  the  modified  UVLM  to  model  LEV  formation  on  finite 
wings.  The  prediction  of  LEV  structure  and  aerodynamic  loads  derived  by  new  modified  UVLM 
are  verified  with  the  observation  of  CFD. 


9.1  Development  of  LEV  structure 

In  this  section,  distribution  of  LEV  sheet  predicted  by  new  LEV  model  is  reported  and  verified 
by  CFD  observation.  As  mentioned  before,  the  motion  of  vortex  sheets  are  determined  by  local 
velocity  field.  For  this  reason,  the  sheet  distribution  can  be  thought  of  as  a  plane  swept  by  local 
streamlines  and  the  figure  of  vortex  sheets  represents  the  flow  structure.  This  section  discusses  the 
accuracy  of  predicted  flow  structure  by  verifying  with  CFD  observation  and  shows  the  development 
of  inner  vortical  structure  of  LEV. 

Figure  9.1  shows  chronological  snapshots  of  LEV  development  for  case  1  (baseline).  Because  the 
flow  is  symmetric,  all  pictures  show  the  left  side  of  wing  and  its  flow  structure.  The  left  and  the 
right  column  of  figure  9.1  respectively  shows  the  observation  of  CFD  and  the  prediction  of  UVLM. 
Herein,  CFD  observation  is  rendered  by  iso-surface  of  Q-criterion  to  depict  the  vortex  structure. 
Although  iso-Q  surface  is  not  quantitatively  equivalent  to  the  vortex  sheet  distribution,  it  allows  for 
qualitative  comparison  of  vortex  flow  structure.  Figure  9.1a  shows  the  initiation  of  LEV  formation. 
The  LEV  sheet  starts  shedding  from  the  leading  edge  around  the  wing  root,  which  is  the  onset 
location  of  LEV  formation,  and  it  corresponds  with  the  CFD  observation.  In  figure  9.1b,  the 
front  line  of  LEV  progresses  to  the  spanwise  direction.  The  predicted  spanwise  progression  of 
LEV  qualitatively  matches  the  CFD  observation  but  UVLM  prediction  is  slightly  delayed.  CFD 
observation  captures  a  few  chordwise  shear  waves  (Kelvin-Helmholtz  wave)  behind  the  leading  edge 
while  UVLM  prediction  shows  a  convection  of  LEV  sheet  with  small  crease  at  the  edge  of  LEV 
sheet,  which  indicates  the  trace  of  the  starting  vortex  of  LEV  shedding.  UVLM  employs  the  cutoff 
function  to  remove  the  singularity  effect  and  small  shear  waves  observed  in  CFD  are  not  simulated 
in  UVLM  prediction.  (Note  that  the  influence  of  shear  waves  are  not  important  compared  to  that 
of  primary  LEV  in  present  work.)  According  to  CFD  observation,  LEV  shedding  region  has  covered 
the  whole  leading  edge  in  figure  9.1c,  however  UVLM  prediction  is  slightly  delayed  and  the  spanwise 
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(b)  t*  =  1.830, 
a  =  28.28  deg. 


(d)  t*  =  2.325, 
a  =  44.53  deg. 


(e)  t*  =  2.505, 
a  =  45.00  deg. 


Figure  9.1:  Comparison  of  LEV  structure  of  CFD  (left)  and  UVLM  (right)  for  case  1. 
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progression  of  LEV  has  not  yet  reached  the  wing  tip.  This  delay  of  spanwise  development  of  LEV 
structure  probably  results  from  the  lack  of  tip  vortices,  which  could  pull  the  surface  flow  structure 
to  the  tip  side.  In  figure  9. Id,  both  CFD  observation  and  UVLM  prediction  shows  the  primary  LEV 
build  up  around  the  leading  edge.  At  the  trailing  edge,  CFD  observation  shows  relatively  large 
shear  wave  and  it  is  also  predicted  by  UVLM.  In  figure  9.1e,  both  CFD  observation  and  UVLM 
prediction  shows  one  developed  primary  LEV  and  several  smaller  waves  on  wing  surface.  CFD 
observation  in  figure  9.1  shows  that  the  tip  vortex  shed  from  the  wing  tip  and  it  directly  connects 
with  the  trailing  vortex  of  TEV.  They  eventually  forms  one  large  and  closed  circle,  whereas  UVLM 
prediction  does  not  have  this  feature  because  UVLM  does  not  simulate  a  detached  or  rolled  up 
tip  vortex.  The  problem  with  simulation  of  the  tip  vortices  in  UVLM  is  discussed  in  subsequent 
section.  Overall,  figure  9.1  provides  a  good  demonstration  of  the  prediction  performance  of  UVLM 
in  the  sense  of  LEV  onset  time,  onset  location,  chordwise  convection,  and  the  basic  pattern  of  LEV 
structure  composed  by  the  primary  LEV.  On  the  other  hand,  it  shows  a  slight  delay  of  spanwise 
development  of  LEV  in  an  early  phase  of  LEV  formation. 

Figure  9.2  shows  the  verification  of  predicted  and  observed  LEV  structure  for  case  6  (tip  twist 
wing).  To  sum  up,  it  can  be  said  that  the  accuracy  of  UVLM  prediction  is  good  to  capture 
the  time-dependent  structural  feature  of  LEV.  As  mentioned  in  previous  section,  one  outstanding 
character  of  twisted  wing  is  that  the  onset  of  LEV  formation  initiates  at  the  semi-half-span  in 
both  left  and  right  wing  and  spanwise  LESP  distribution  has  an  M-shaped  profile.  As  shown  in 
figure  9.2a,  LEV  onset  mark  is  indicated  at  around  2 y/b  =  0.6  in  both  CFD  observation  and  UVLM 
prediction.  In  figure  9.2b,  the  span  of  LEV  shedding  area  rapidly  progresses  towards  wing  root  and 
wing  tip  and  covers  the  whole  leading  edge  in  CFD  observation.  On  one  hand,  the  UVLM  prediction 
shows  that  spanwise  progression  of  the  LEV  has  reached  the  wing  root,  on  the  other  hand,  it  delays 
toward  the  wing  tip,  similar  to  the  verification  of  case  1.  As  discussed  in  previous  section,  the 
onset  of  LEV  formation  of  the  twisted  wing  starts  from  the  semi- half- span.  The  strength  of  the 
leading  edge  suction  has  twin  peak  (called  M-shaped  profile)  at  the  point  and  the  trace  of  profile 
is  preserved  for  a  few  seconds.  As  a  result,  the  primary  LEV  forms  a  slightly  distorted  shape, 
which  would  be  illustrated  as  a  hyperbolic  hyperboloid  in  full  wing  model,  as  shown  in  figure  9.2e. 
(This  hyperboloid  structure  is  easily  understood  by  comparing  the  line  of  the  leading  edge  and  the 
outline  of  primary  LEV.  )  UVLM  correctly  predicts  the  structural  feature  of  LEV  for  the  twisted 
wing.  The  M-shaped  LEV  structure  seems  to  be  transient  and  soon  converges  into  one  cylindrical 
LEV  soon  after  a  second.  According  to  UVLM  prediction  at  t*  =  3.300  as  shown  in  figure  9.3,  it 
can  be  confirmed  that  the  hyperboloid  structure  has  been  absorbed  by  new  major  cylindrical  LEV 
and  it  finally  attains  the  same  LEV  structure  as  that  of  baseline  case. 

Figure  9.4  shows  the  comparison  of  LEV  structure  of  CFD  observation  and  UVLM  prediction 
through  the  time  sequence.  As  mentioned  in  previous  section,  the  onset  of  LEV  in  the  swept 
wing  occurs  near  the  wing  tip.  This  behavior  is  confirmed  through  figures  9.4a  and  9.4b.  UVLM 
prediction  captures  well  with  CFD  observation  for  the  LEV  formation  and  development  towards 
the  wing  root,  as  seen  in  figures  9.4b  and  9.4c.  Finally,  as  shown  in  figures  9.4d  and  9.4e,  the 
LEV  forms  a  conical  shape  which  is  well  known  in  delta  wings.  To  sum  up,  UVLM  successfully 
predicts  the  LEV  structure  compared  to  CFD  observation  in  case  10.  However,  CFD  also  shows 
the  interaction  between  LEV  and  the  tip  vortices.  Figure  9.5  is  the  same  snapshot  of  figure  9.4e 
rendered  from  the  front-left  viewpoint.  According  to  figure  9.5,  primary  LEV  shows  a  conical  shape 
from  the  wing  root  (left  side  of  the  picture)  to  the  wing  tip  (bottom-right  corner  of  the  picture)  and 
contacts  the  tip  vortex.  Considering  Helmholtz  law,  LEV  would  feed  itself  to  the  tip  vortex  there 
and  they  eventually  connect  with  each  other  at  the  corner  of  the  wing  tip.  It  could  be  geometrically 
important  feature  to  simulate  in  UVLM  prediction  and  is  discussed  in  subsequent  section. 
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(b)  t*  =  1.710, 
a  =  24.15  deg. 


(c)  t*  =  2.010, 
a  =  34.46  deg. 


(d)  t*  =  2.280, 
a  =  43.55  deg. 


(e)  t*  =  2.505, 
a  =  45.00  deg. 


Figure  9.2:  Comparison  of  LEV  structure  of  CFD  (left)  and  UVLM  (right)  for  case  6. 
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Lastly,  the  result  of  UVLM  prediction  of  LEV  development  with  inner  structure  is  presented.  Fig¬ 
ure  9.6  shows  a  series  of  half-cut  snapshots  of  LEV  sheet  distribution  predicted  by  UVLM  for  case  1 
(baseline).  UVLM  prediction  successfully  captures  the  process  of  LEV  formation;  the  onset  of  LEV 
formation  initiates  from  the  wing  root  (figure  9.6a),  the  LEV  shedding  area  progresses  to  span- 
wise  direction  and  the  shed  LEV  sheet  convects  to  downstream  (figure  9.6b),  LEV  structure  starts 
building  up  (figure  9.6c),  an  reattached  region  of  LEV  sheet  retreats  (moves  forward)  at  around 
mid-chord  (figure  9.6d),  a  first  inner  turn  is  formed  in  the  primary  LEV  structure  (figure  9.6e), 
and  major  vortical  structure  is  finally  developed  on  the  wing  upper  surface  (figure  9.6f).  Though 
a  verification  of  inner  structure  of  vortex  is  technically  difficult  in  present  work,  the  prediction 
result  of  UVLM  seems  reasonable,  considering  several  fundamental  laws  and  some  experimental 
observation  of  LEV  structure,  for  example  the  reference.57 


9.2  Prediction  of  aerodynamic  loads  with  LEV 

Verification  of  lift  coefficient  Cl  and  momentum  coefficient  about  aerodynamic  center  Cm  for  case 
1  (baseline)  is  presented  in  figure  9.7.  When  the  wing  starts  pitching  at  around  t*  =  1,  Cl  increases 
sharply  due  to  the  effect  of  apparent  mass  and  Cm  drops  because  of  the  moment  of  apparent  mass  is 
loaded  about  the  half  chord.  Predicted  Cm  is  overestimated  compared  to  CFD  observation  during 
this  period  because  UVLM  is  based  on  inviscid  flow  model.  The  apparent  mass  subsides  after  the 
end  of  pitch  acceleration  and  Cl  rises  as  the  angle  of  attack  a  increases.  Cm  of  CFD  observation 
gradually  decreases  with  increase  of  a  while  the  Cm  in  UVLM  prediction  is  nearly  constant  because 
UVLM  is  based  on  potential  flow  theory  and  Cm  is  defined  at  aerodynamic  center.  After  the  onset 
of  LEV  formation,  both  observed  Cl  and  predicted  Cl  keeps  increasing  but  the  rate  of  increase 
gets  weaker  in  UVLM  prediction  whereas  the  slope  of  Cl  of  CFD  observation  does  not  change 
through  the  LEV  initiation.  This  difference  in  Cl  slope  eventually  results  in  a  maximum  of  about 
15%  discrepancy  in  Cl -  For  Cm,  the  predicted  Cm  turns  its  variation  from  up  to  down  at  the  onset 
of  LEV  formation  and  the  trend  becomes  the  same  as  CFD  observation  after  LEV  formation.  At 
around  t*  =  2.3,  the  pitching  wing  decelerates  and  both  predicted  and  observed  Cl  drop  due  to 
negative  apparent  mass.  On  the  other  hand,  due  to  the  location  of  pitch  axis,  this  negative  apparent 
mass  raises  the  Cm-  Since  the  wing  attains  the  end  of  motion  (a  =  45  deg)  and  stops  pitching, 
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(a)  t*  =  1.080, 
a  =  2.54  deg. 


(b)  t*  =  1.275, 
a  =  9.20  deg. 


(c)  t*  =  1.620, 
a  =  21.06  deg. 


(d)  t*  =  1.800, 
a  =  27.24  deg. 


(e)  t*  =  2.145, 
a  =  39.10  deg. 


Figure  9.4:  Comparison  of  LEV  structure  of  CFD  (left)  and  UVLM  (right)  for  case  10. 
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Figure  9.5:  Interaction  of  LEV  and  tip  vortex,  t*  =  2.145,  a  =  39.10  deg. 


Cl  decreases  by  development  of  LEV  on  wing  surface.  According  to  figure  9.6,  LEV  keeps  growing 
and  moves  downstream.  Thus,  Cm  in  UVLM  prediction  decreases  after  the  wing  stops  pitching.  A 
fluctuation  is  found  in  predicted  Cl  and  Cm  for  this  period  and  it  can  be  thought  as  an  influence  of 
instability  waves  on  the  wing  surface.  As  the  primary  LEV  develops  its  structure,  the  lower  surface 
of  LEV  pushes  the  circulation  flow  to  the  wing  surface.  The  circulation  has  Kelvin- Hehxdioltz  wave 
due  to  shear  force  near  the  surface  and  they  are  observed  as  a  fluctuation  pattern  in  Cl  and  Cm- 

Analysis  of  the  aforementioned  discrepancy  in  Cl  about  LEV  formation  has  centerd  on  three 
probable  causes.  First  probable  cause  is  the  underestimation  of  the  leading  edge  suction  force  Fs. 
LESP  concept  postulates  that  the  flow  around  the  leading  edge  must  be  preserved  under  LEV 
formation.  Then,  the  value  of  LESP(t)  plateaus  at  LESPcrit  during  LEV  shedding.  Figure  9.8 
is  a  comparison  of  maximum  LESP  between  UVLM  and  CFD.  The  method  of  conversion  from 
CFD  data  to  LESP  is  introduced  by  Narsipur  et  al .126  In  figure  9.8,  CFD  observation  shows  that 
the  leading  edge  suction  does  not  necessarily  plateau  since  LEV  initiation  and  keeps  increasing  for 
a  while,  whereas  LESPmax(f)  of  UVLM  forced  to  be  plateaued.  The  period  of  this  discrepancy 
in  LESP  between  UVLM  prediction  and  CFD  observation  corresponds  with  the  period  of  high 
incidence  (a  ~  45  deg)  in  which  the  suction  force  Fs  becomes  predominant  in  Cl- 

Second  probable  cause  is  the  fact  that  current  UVLM  does  not  simulate  the  influence  of  tip  vortices. 
Current  UVLM  does  have  a  tip  vortex  model  as  an  option  but  it  is  not  usually  activated  because 
it  is  still  under  development.  In  addition,  previous  case  studies  showed  that  tip  vortices  did  not 
largely  affect  the  prediction  of  LEV  onset  time  and  location  except  for  a  short  span  wing.  However, 
tip  vortices  would  positively  influence  the  forces  on  the  lifting  surface  to  some  extent,  even  if  the 
wing  geometry  has  broad  wing  span.  Thus,  the  improvement  of  the  tip  vortex  model  could  reduce 
the  discrepancy  in  Cl  after  LEV  formation. 

Third  probable  cause  is  the  delay  of  development  of  LEV.  As  be  seen  in  previous  section,  UVLM 
prediction  shows  qualitatively  good  agreement  with  CFD  observation.  However,  it  is  also  confirmed 
that  the  prediction  of  LEV  development  in  spanwise  direction  is  somewhat  slower  than  the  result 
of  CFD  observation.  The  delay  in  spanwise  progression  of  LEV  onset  must  cause  underestimation 
of  the  influence  of  LEV.  In  order  to  assess  the  LEV  model  quantitatively,  one  more  verification 
is  implemented  about  the  strength  of  circulation  variation,  using  LESP-modulated  discrete-vortex 
method  (LDVM,  see  the  reference2).  In  this  investigation,  because  LDVM  is  two-dimensional 
flow  solution,  current  UVLM  calculates  for  a  pseudo-infinite  wing  (AR  =  20,000)  to  simulate 
approximately  two-dimensional  flow  field.  (This  method  is  introduced  in  the  reference.127) 
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(a)  t*  =  1.710,  a  =  19.51  deg. 


(b)  t*  =  2.040,  a  =  28.28  deg. 


(c)  t*  =  2.280,  a  =  43.55  deg. 


(d)  t*  =  2.565,  a  =  45.00  deg. 


(e)  t*  =  2.820,  a  =  45.00  deg. 


(f)  t*  =  3.300,  a  =  45.00  deg. 


Figure  9.6:  Time  variation  of  development  of  LEV  structure  (case  1). 
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Lift  coefficient  (CL) 


0  0.5  1  1.5  2  2.5  3 

Nondimensional  time  (t*) 


Figure  9.8:  Comparison  of  LESPmax(t )  for  case  1. 
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Angle  of  attack  (Alpha),  deg 


Figure  9.9:  Comparison  of  circulation  strength  for  case  Ex  1. 


The  result  of  comparison  of  UVLM  for  a  pseudo-infinite  wing  and  LDVM  is  presented  in  figure  9.9. 
Following  the  onset  of  LEV  formation,  both  UVLM  and  LDVM  show  that  the  wing  bound  circu¬ 
lation  decreases  while  the  strength  of  LEV  increases.  On  the  other  hand,  LEV  formation  does  not 
clearly  influence  the  increase  rate  of  TEV.  UVLM  shows  the  same  trend  in  these  three  circulations 
but  the  strengths  of  LEV  and  TEV  are  weaker  than  that  of  LDVM  data.  The  difference  of  devel¬ 
opment  of  circulation  must  affect  the  distribution  of  LEV  sheet  in  UVLM.  In  detail,  the  stronger 
LEV  sheet  must  roll  up  faster  and  be  slower  to  convect  downstream.  Consequently,  the  sheet  dis¬ 
tribution  must  be  more  dense  around  the  leading  edge.  Therefore,  the  dense  LEV  sheet  interacts 
with  stronger  circulation  on  the  lifting  surface.  The  concentrated  LEV  structure  developed  by  fast 
LEV  growth  would  mitigate  the  deficiency  of  wing  bound  circulation  by  LEV  formation.  For  this 
reason,  the  weak  strength  of  LEV  in  UVLM  prediction  would  cause  the  discrepancy  in  Cl- 

These  three  probable  causes  do  not  completely  explain  all  the  symptoms.  They  also  relate  each 
other.  To  conclude,  it  should  be  thought  that  the  further  improvement  of  accuracy  of  prediction  of 
aerodynamic  loads  requires  that  these  factors  be  addressed. 


9.3  Computation  time 


Table  9.1:  Computation  time  for  baseline  case. 


Code 

Computation  time 

UVLM  with  LEV  model 

~  25  minutes 

REACTMB 

~  1  day 

Result  of  benchmark  test  of  computation  time  is  summarized  in  table  9.1.  UVLM  is  a  representative 
low-order  method  based  on  theoretical  model  and  NCSU’s  REACTMB  is  taken  as  a  referential 
cutting-edge  CFD  method  based  on  incompressible  Navier-Stokes  solver.  The  benchmark  for  both 
the  codes  is  measured  by  the  same  computer  environment  (NC  State  HPC  as  of  2015  -  2016).  The 
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(a)  Prior  to  amalgamation  (t*  =  3.045).  (b)  After  amalgamation  (t*  =  3.060). 

Figure  9.10:  Snapshots  of  three-dimensional  amalgamation  method 


benchmark  test  shows  that  current  UVLM  run  about  50  times  as  faster  as  REACTMB. 

Table  9.2:  Contribution  of  LEV  model  and  vortex  roll  up  to  computation  time. 


Test  No. 

Combination 

LEV  model 

Roll  up 

cpu  time 

1 

UVLM  without  roll  up 

off 

off 

24  sec 

2 

Conventional  UVLM 

off 

on 

93  sec 

3 

UVLM  with  LEV  model 

on 

on 

168  sec 

Conventional  UVLM  is  neither  postulated  to  be  used  for  simulating  quick  motion  kinematics  nor  in 
high  angle  of  attack  condition.76  Conventional  UVLM,  of  course,  does  not  simulate  LEV.  Therefore, 
the  influence  of  including  LEV  model  to  the  total  computation  time  is  of  interest.  The  table  9.2 
shows  a  summary  of  problem  analysis  for  reduction  of  computation  time,  which  compares  the 
computation  time  (cpu  time)  for  a  few  combinations  of  LEV  model  and  vortex  roll  up.  This 
benchmark  was  run  the  same  condition  of  case  1  using  low  resolution  discretization  (TV  0(10)). 

The  most  important  information  for  present  work  is  the  increase  of  computation  time  by  addition 
of  the  LEV  model.  Comparing  the  combination  test  2  to  3,  it  is  known  that  the  new  LEV  model 
increases  +80.6%  of  computation  time.  However,  this  is  not  mainly  caused  by  LEV  model  itself 
but  by  vortex  roll  up.  Taking  combination  test  1  and  2  in  table  9.2  as  an  example,  it  is  shown  that 
vortex  roll  up  occupies  74%  (93-24=69  sec)  of  conventional  UVLM  run  (93sec).  This  is  because 
vortex  roll  up  needs  to  estimate  influence  coefficients  for  all  vortex  rings  in  every  time- step.  The 
computational  time  of  vortex  roll  up  is  proportional  to  the  square  of  total  number  of  vortex  rings 
and  the  LEV  model  eventually  adds  LEV  rings  in  computation.  Thus,  an  engineering  item  to 
reduce  the  number  of  vortex  elements  would  be  required  to  save  the  computation  time.  As  pointed 
out  by  Suresh  Babu  et  a/.,106  the  most  promissing  avenue  for  reducing  the  number  of  LEV  grids  is 
an  amalgamation  of  shed  vortex  elements. 


109 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


9.4  Demonstration  of  three-dimensional  amalgamation  method 

As  mentioned  in  previous  section,  main  objection  of  amalgamation  in  vortex  method  is  to  reduce 
the  computation  time  and  to  prevent  impingement  of  the  vortex  sheets.  Current  UVLM  has  an 
amalgamation  model  for  LEV  sheet  and  the  test  result  is  presented  in  figure  9.10.  The  present 
amalgamation  method  simply  amalgamates  the  inner  turn  of  primary  LEV  by  merging  the  spanwise 
row  of  LEV  sheets  when  the  external  angle  of  LEV  surface  projected  on  XZ  plane  (y  =  0)  0 
exceeds  the  threshold  6a  =  180  deg.  (For  convenience  of  presentation,  the  amalgamation  condition 
is  turned  “off”  for  a  while,  to  allow  some  build  up  of  the  LEV,  and  then  turned  “on”  to  merge 
the  accumulated  inner  turns.  )  As  be  seen,  an  inner  turn  structure  of  primary  LEV  shown  in 
figure  9.10a  is  amalgamated  into  a  single  line  in  figure  9.10b  and  the  outer  half  circular  surface  of 
LEV  sheet  is  left  because  0a  =  180.  The  present  action  of  amalgamation  has  collapsed  28  rows  of 
LEV  sheet  in  the  primary  LEV  and  eventually  1680  LEV  panels  has  been  eliminated  between  the 
time-step.  As  a  result,  correct  operation  of  amalgamation  was  confirmed  in  the  test  and  it  can  be 
said  the  basic  idea  of  amalgamation  method  for  vortex  ring  discretization  works  well.  Unfortunately 
clear  contribution  to  reduction  of  computation  time  was  not  shown  in  this  test  condition  because 
the  amalgamation  was  always  activated  toward  end  of  the  simulation  sequence.  However  a  large 
saving  would  be  expected  in  general  situation.  On  the  other  hand,  current  amalgamation  method 
indicates  a  side  effect  and  it  is  discussed  in  following  section. 


9.5  Known  problems 


While  current  UVLM  has  demonstrated  adequate  prediction  of  LEV  formation,  a  few  numeric 
problems  have  also  been  found. 

9.5.1  Geometric  restrictions 

UVLM  deals  with  a  hyperplane  problem,  and  three-dimensional  flow  field  is  determined  by  lumped 
vortex  sheets  which  form  a  curvilinear  surface.  Therefore,  UVLM  is  governed  by  not  only  fluid 
dynamics  but  also  geometry.  The  problem  is  that  the  vortex  structure  which  determines  the  flow 
field  is  not  always  available  from  vortex  sheet  distribution.  UVLM  may  not  be  acceptable  in 
particular  situations  which  violate  geometric  restrictions. 

First  problem  of  geometric  restriction  in  UVLM  is  the  shape  of  vortex  panel  element.  Vortex  method 
commonly  defines  a  unit  vortex  element  as  quadrilateral  form.  Quadrilateral  has  four  corners  and 
can  bend  in  one  axis  including  two  of  four  corners.  This  geometric  character  of  quadrilateral  vortex 
panel  does  not  allow  bends  in  two  direction  if  it  is  in  a  cross  point  of  two  rolled  structures.  For 
example,  UVLM  cannot  correctly  calculate  the  vortex  structure  at  the  region  where  conical  LEV 
turns  and  trailing  tip  vortex,  as  be  seen  in  case  10.  Quadrilateral  vortex  element  is  not  always 
suitable  for  describing  a  complicated  flow  structure. 

Second  problem  of  geometric  restriction  is  about  resolution.  To  form  a  curvilinear  surface  by  two 
or  less  points  is  impossible.  Similarly,  the  vortical  structure  requires  to  have  sufficient  number  of 
vortex  sheets  to  be  described  correctly.  However,  the  number  of  vortex  sheets  depends  on  the  flow 
field  and  numeric  condition  such  as  time-step  and  panel  size.  Furthermore,  LEV  and  tip  vortices 
is  known  to  have  a  multiple  vortex  structure  but  the  limit  of  resolution  cannot  always  depict 
them.  The  formation  of  secondary  and  subsequent  vortices  often  causes  an  unreal  distortion  and 
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impingement  in  vortex  sheet.  Thus,  regionally  low  resolution  is  often  seen,  resulting  in  difficulty  in 
simulating  correct  vortex  flow  structure,  and  causes  a  numeric  error  in  UVLM.31 

Third  problem  of  geometric  restriction  is  provided  by  singularity  of  vortex  element  itself.  It  is 
known  that  shed  vortex  grids  gradually  congest  and  form  a  vortical  structure.105  Accordingly,  they 
automatically  approach  to  their  singularity  zone.  As  mentioned  before,  application  of  cutoff  radius 
can  mitigate  an  individual  input  error  but  cannot  prevent  from  the  tendency  of  congestion  which 
would  increase  the  input  error.  Singularity  is  also  problematic  in  a  situation  which  does  not  have 
enough  space  to  build  a  vortex  structure  by  discretized  vortex  elements.  For  one  example,  a  sharp 
leading  edge  wing  generates  LEV  in  low  incidence  angle  but  this  low  incidence  situation  does  not 
give  a  sufficient  space  of  recirculation  region  behind  LEV  to  develop.  For  another  example,  a  case 
which  is  very  low  LESPcrit  and  slow  K  motion  tends  to  be  unstable  because  LEV  sheds  in  a  narrow 
space  between  bulk  flow  and  the  wing  surface  and  easily  enter  the  singularity  zone.  Thus,  time-step 
and  cutoff  radius  also  restrict  the  case. 

Essentially,  these  geometric  restrictions  are  mathematically  fundamental  problem  of  conventional 
vortex  roll  up.  Process  of  vortex  roll  up  is  a  continuous  mapping  from  last  time-step  to  present. 
The  problem  is  that  the  ordinary  mapping  function  in  vortex  method  must  include  singularities  of 
vortex  filaments.  This  is  one  of  difficulties  of  all  vortex  methods.  Consider  a  sphere  covered  by  a 
number  of  quadrilateral  elements  as  an  example.  In  an  ordinary  mapping,  it  is  possible  to  fit  each 
quadrilateral  element  for  the  curvature  and  the  quadrilaterals  at  the  poles  must  be  deformed  to  al¬ 
most  triangle  shape.  Singularity  restriction  in  vortex  method  could  not  allow  this  mapping.  Several 
experimental  observation  shows  that  LEV  shed  from  a  finite  leading  edge  gradually  transforms  to 
be  a  curvilinear  shape.  However,  current  LEV  model  in  UVLM  does  not  always  provide  sufficient 
degree  of  freedom  to  form  an  arbitrary  vortex  structure  because  of  these  geometric  restrictions  as 
already  seen.  Furthermore,  the  actual  vortex  structure  does  not  always  satisfy  the  mathemati¬ 
cal  consistency.  For  example,  if  the  problem  requires  to  simulate  a  long-term  wake  distribution, 
UVLM  must  be  able  to  model  vortex  dissipation,  separation  and  reconnection,  as  introduced  in 
the  reference.128  Vortex  reconnection  is  an  intriguing  topic  because  it  is  a  violation  of  topological 
fluid  dynamics.  This  phenomenon  is  thought  to  be  caused  by  viscous  interaction  and  turbulent 
transition  but  UVLM  tends  to  maintain  the  topological  consistency  and  preserve  the  helicity  in  the 
system  because  it  is  based  on  potential  flow  theory. 

One  of  solution  of  geometric  difficulty  is  an  amalgamation  method  as  shown  in  previous  section. 
However,  it  is  also  noticed  that  the  amalgamation  method  is  fundamentally  nonlinear  operation 
and  violation  of  topology,  while  amalgamation  method  has  to  be  identified  to  resulting  velocity 
field.  To  prevent  the  topological  violation  and  ensure  the  identification  of  velocity  field,  UVLM 
must  have  sufficient  degree  of  freedom  but  it  is  disturbed  by  aforementioned  geometric  restrictions. 
Amalgamation  model  in  current  UVLM  does  not  have  clear  solution  to  cover  all  these  problems  and 
becomes  unstable  after  the  operation  of  amalgamation.  To  sum  up,  further  improvement  of  LEV 
model  will  be  required  to  simulate  more  advanced  conditions  and  the  solution  must  be  a  technique 
to  satisfy  or  avoid  these  geometric  difficulties. 

9.5.2  Management  of  numeric  instability 

Generally,  control  of  the  UVLM  numeric  instability  is  difficult  to  manage  due  to  nonlinearity  of 
modelling.  Thus,  the  solution  of  instability  in  UVLM  has  tended  to  be  ad  hoc.  Current  UVLM 
and  LEV  model  does  not  equip  any  control  method  against  the  instability  and  there  are  no  ways 
to  mitigate  the  amplification  of  numeric  error. 
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The  formulation  of  velocity  field  of  UVLM  is  not  suitable  for  a  common  stability  analysis.  Procedure 
of  the  common  stability  analysis  consists  of  to  give  a  perturbation  (in  most  cases,  the  perturbation 
is  defined  as  eat+lkx  where  a  is  an  amplifier,  i  is  imaginary  unit,  k  is  a  wave  number),  to  reduce 
the  original  equation  to  an  eigenfunction  by  linealization,  to  evaluate  the  distribution  of  the  eigen¬ 
values  (amplifier).  One  example  of  instability  analysis  in  vortex  method  was  reported  by  Crow.128 
Crow  modelled  the  two  trailing  vortices  from  wing  tip  as  a  couple  of  parallel  vortex  filaments  and 
implemented  an  eigen- mode  analysis.  In  his  instability  analysis,  the  two  tip  vortices  were  assumed 
to  locate  in  sufficient  distance  from  each  other  and  both  tip  vortex  can  be  linearized  as  an  cylindri¬ 
cal  velocity  field  formulated  by  Bessel  function.  However,  the  nonlinearity  of  UVLM  system  does 
not  allow  implementation  of  linealization  because  the  relative  distance  of  discrete  vortices  varies  in 
broad  range. 

There  were  a  number  of  instability  analysis  for  UVLM  from  a  different  perspective.  Subvortex 
technique  proposed  by  Maskew129  was  an  example  of  uncommon  instability  analysis  for  UVLM.  In 
subvortex  technique,  the  numeric  instability  was  simply  evaluated  by  a  function  of  the  distance  of 
vortex  filaments.  This  function  shows  an  optimal  vortex  panel  size  and  the  unstable  vortex  panels 
were  temporally  split  by  a  few  subvortices  to  avoid  the  numeric  error.  Subvortex  technique  equips 
the  method  of  instability  analysis  in  relative  location  but  the  derivation  of  the  optimal  subvortex 
size  is  cumbersome  and  not  suitable  for  three-dimensional  UVLM. 

To  sum  up,  the  tendency  of  the  instability  is  not  clearly  determined  in  UVLM.  With  a  new  idea 
of  stability  analysis  method,  some  numeric  errors  in  UVLM  would  potentially  be  predicted  and 
appropriately  inhibited  in  future. 
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Chapter  10 


Conclusions  and  Future  Work 


The  major  focus  of  the  current  work  was  to  extend  the  LESP  concept,  which  provides  a  measure  of 
leading-edge  suction  in  two-dimensional  geometries,  to  handle  low-order  prediction  of  LEV  initiation 
and  formation  on  three-dimensional  finite- wing  systems.  Earlier  work  on  airfoils  in  2D  unsteady 
flow  showed  that  when  the  instantaneous  LESP  exceeds  a  critical  value,  LEV  shedding  occurs. 
One  of  the  major  insights  from  that  research  effort  was  that  the  critical  LESP  is  dependent  only 
on  the  airfoil  and  Reynolds  number,  and  is  largely  independent  of  motion  kinematics  so  long  as 
the  LEV  formation  is  not  preceded  by  significant  trailing-edge  separation.  Thus,  for  this  class  of 
motions,  the  critical  LESP  for  a  given  airfoil  and  Reynolds  number  can  be  determined  from  CFD 
or  experiment  for  one  prototypical  motion,  and  can  then  be  used  for  any  other  high-rate  motion, 
including  arbitrary  combinations  of  pitch,  plunge,  and  surge. 

The  current  work  started  with  the  hypothesis  that  the  flow  physics  governing  LEV  intiation  on 
airfoil  and  finite- wing  geometries  are  the  same,  provided  that  the  induced- flow  effects  due  to  finite- 
wing  effects  are  taken  into  consideration.  The  implication  is  that  the  value  of  critical  LESP  for  a 
given  airfoil  and  Reynolds  number  will  also  be  applicable  to  a  finite  wing  which  uses  that  airfoil 
geometry.  To  test  this  hypothesis  about  LEV  initiation  on  finite  wings,  a  standard  unsteady  vortex 
lattice  method  (UVLM)  was  adapted  to  calculate  the  spanwise  distribution  of  instantaneous  LESP 
along  the  wing  at  every  time  instant  in  the  motion.  LEV  initiation,  from  low-order  prediction,  is 
assumed  to  occur  at  the  time  instant  and  spanwise  location  when  the  local  value  of  LESP  equals 
the  two-dimensional  value  of  critical  LESP.  Low-order  predictions  for  LEV  initiation  for  several 
wing  planforms  (13  cases  with  variations  in  planform  geometry,  pivot  location,  pitch  rate,  twist 
angle,  aspect  ratio,  and  cross-section  shape)  were  compared  with  RANS  CFD  predictions. 

It  was  seen  that,  for  any  given  airfoil  shape  and  Reynolds  number,  the  low-order  predictions  of  angle 
of  attack  and  spanwise  location  of  LEV  initiation  were  acceptably  close  to  the  corresponding  CFD 
predictions.  Specifically,  the  angle  of  attack  for  LEV  initiation  between  the  two  approaches  were 
within  zb  2-degree  accuracy  against  a  40-degree  spread  in  angle  of  attack  for  all  the  cases  considered. 
The  current  results  confirm  that  initiation  of  LEV  formation  on  finite  wings  is  governed  by  the 
same  critical  value  of  the  LESP  that  is  applicable  in  2D  flow.  The  major  implication  is  that  the 
critical  value  of  LESP  obtained  from  a  two-dimensional  study  for  one  “high-rate”  motion  can  be 
used  for  prediction  of  LEV  initiation  on  any  finite- wing  geometry  undergoing  any  other  “high-rate” 
motion.  A  further  implication  is  that  the  LESP  concept  opens  up  avenues  for  LEV  flow  control  by 
tailoring  either  the  instantaneous  or  critical  LESP  at  different  portions  of  the  wing  span. 

The  second  part  of  the  research  focused  on  modifying  the  UVLM  to  handle  LEV  formation.  In  this 
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modified  UVLM,  a  vortex  sheet  from  the  leading  edge  is  modeled,  and  its  geometry  is  calculated  in 
each  time  step  by  assuming  that  the  corner  points  of  the  panels  forming  the  sheet  convect  with  the 
local  velocity.  An  important  element  of  this  effort  is  the  development  of  a  new  idea — pseudo  vortex 
ring — for  vortex  lattice  method  to  enable  calculation  of  LEV  formation  in  spite  of  problems  with 
vortex-ring  discretization  at  the  leading  edge.  The  modified  UVLM  was  successfully  in  achieving 
good  prediction  of  time-dependent  development  of  LEV  structure.  For  the  calculation  of  aerody¬ 
namic  loads,  an  approximately  15%  of  underestimation  of  circulation  strength  of  LEV  was  observed 
in  current  low-order  LEV  formation  model.  On  the  other  hand,  the  low-order  computations  are 
50  times  faster  than  the  RANS  CFD  computations,  which  makes  the  current  method  suitable  for 
engineering  and  design  tasks. 

While  the  low-order  prediction  of  LEV  formation  shows  promise,  significant  challenges  remain  in 
modeling  of  the  sheet  roll  up  and  its  intersection  with  the  wing  geometry.  Discrete- vortex  amalga¬ 
mation  techniques,  studied  in  the  current  effort  for  2D  flows,  shows  some  promise  for  reducing  the 
number  of  vortex  lattices  in  LEV  sheets  from  finite  wings.  It  is  likely  that  some  adaption  of  vortex 
amalgamation  applied  to  vortex  sheets  from  leading  edges  will  allow  for  robust  modeling  of  LEV 
shedding  on  finite  wings.  Improved  modeling  of  the  rolled- up  tip  vortex  structures  is  also  necessary 
if  prediction  accuracy  for  low  aspect-ratio-wings  is  to  be  improved. 

Suggested  future  work  includes  (i)  augmentation  of  the  current  UVLM  formulation  to  include  better 
modeling  of  tip-vortex  shedding,  (ii)  development  of  a  hybrid  vortex-sheet  /  vortex- filament  model 
for  LEV  shedding  from  finite  wings,  and  (iii)  extensions  to  the  UVLM  to  handle  LEV  formation, 
growth,  and  detachment,  allowing  for  modeling  of  intermittent  LEV  shedding. 
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Abstract 

The  aim  of  this  research  effort  was  to  extend  earlier  work  on  two-dimensional  airfoils  undergoing  high- 
intensity  unsteady  motions  dominated  by  leading-edge  vortex  (LEV)  shedding  to  finite-wing  geometries.  In 
the  earlier  work  (supported  by  AFOSR  grant  FA9550-1 0-1  -0120,  Pis:  Gopalarathnam,  Edwards,  and  Ol, 

PM:  Dr.  Douglas  Smith),  an  integrated  theoretical,  computational,  and  experimental  research  effort  was 
used,  in  which  experiments  and  higher-order  computations  were  used  to  develop  a  low-order  method  for 
unsteady  aerodynamic  analysis  of  airfoils  and  plates  undergoing  large-amplitude,  high-rate  motions. 

Owing  to  the  fact  that  such  flows  are  dominated  by  LEV  shedding  and  flow  separation,  they  are  well  outside 
the  validity  of  classical  theoretical  methods. 

The  major  contribution  from  the  earlier  research  effort  was  in  identifying  the  importance  of  leading-edge 
suction  in  governing  the  initiation,  formation,  and  termination  of  vortex  shedding  from  rounded  leading 
edges  of  unsteady  airfoils  in  two-dimensional  flow.  It  was  shown  that  the  value  of  this  leading-edge  suction 
at  any  time  instant  in  unsteady  flow  could  be  tracked  in  an  unsteady  thin  airfoil  theory  using  an  inviscid 

parameter  which  was  named  the  Leading-Edge  Suction  Parameter,  or  LESP.  When  the  instantaneous 
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LESP  exceeds  a  critical  value,  LEV  shedding  occurs.  One  of  the  major  insights  from  that  research  effort 
was  that  the  critical  LESP  is  dependent  only  on  the  airfoil  and  Reynolds  number,  and  is  largely 
independent  of  motion  kinematics  so  long  as  the  LEV  formation  is  not  preceded  by  significant  trailing-edge 
separation.  Typically,  LEV  formation  without  accompanying  trailing-edge  separation  occurs  in  high-rate 
and  high-reduced-frequency  motions.  Thus,  for  this  class  of  motions,  the  critical  LESP  for  a  given  airfoil  and 
Reynolds  number  can  be  determined  from  CFD  or  experiment  for  one  prototypical  motion,  and  can  then  be 
used  for  any  other  high-rate  motion,  including  arbitrary  combinations  of  pitch,  plunge,  and  surge.  This 
insight  was  used  to  augment  an  inviscid  unsteady  thin  airfoil  theory  (Ramesh  et  al.,  Theor.  Comput.  Fluid 
Dyn.,  January  2013,  DO1 1 0.1 007/s001 62-012-0292-8)  with  the  addition  of  a  discrete-vortex  method  to 
handle  intermittent  LEV  shedding.  The  resulting  method,  named  the  LESP-modulated  discrete  vortex 
method,  or  LDVM,  is  described  in  detail  in  Ramesh  et  al.  (Journal  of  Fluid  Mechanics,  Volume  751 ,  July 
2014,  pp  500-538.  DOI:  http://dx.doi.Org/1 0.1 01 7/jfm.201 4.297).  The  results  from  LDVM  show  remarkably 
good  agreement  in  forces  and  flow  fields  with  computational  fluid  dynamics  (CFD)  and  experiments.  The 
advantage  of  the  rapid  computational  capability  offered  by  a  low-order  method  like  the  LDVM  was 
illustrated  in  the  prediction  of  non-linear  aeroelastic  behavior  of  airfoils  having  LEV  shedding.  As  described 
in  Ramesh  et  al.  ( Journal  of  Fluids  and  Structures,  Vol.  55,  May  201 5,  pp.  84-1 05. 
http://dx.doi.Org/1 0.1 01 6/j.jfluidstructs.201 5.02.005),  such  aeroelastic  analysis  requires  computation  for  at 
least  1 ,000~convective  times  to  establish  the  long-time  aeroelastic  behavior,  which  would  be  prohibitive 
with  any  high-order  CFD-like  method. 

Motivated  by  the  success  of  the  two-dimensional-flow  work,  the  research  in  the  current  effort  was  aimed  at 
extending  a  three-dimensional  inviscid  analysis  method  to  handle  vortex  shedding  from  rounded  leading 
edges  of  finite  wings.  Systematic  studies  of  finite  wings  with  LEV  shedding  were  analyzed  using  an 
incompressible  Reynolds-Averaged  Navier  Stokes  (RANS)  CFD  method,  the  results  of  which  were  used  for 
the  development  of  the  low-order  method.  At  the  foundation  of  the  low-order  finite-wing  analysis  method  is 
a  traditional  unsteady  vortex  lattice  method  (UVLM).  In  the  first  part  of  the  current  work,  the  objective  was  to 
assess  the  applicability  of  the  LESP  concept  to  prediction  of  the  time  instant  and  spanwise  location  of  LEV 
initiation  on  a  finite  wing  undergoing  unsteady  motion.  The  UVLM  was  used  to  calculate  the  spanwise 
variation  of  LESP  at  every  time  step.  LEV  initiation,  from  low-order  prediction,  is  assumed  to  occur  at  the 
time  instant  and  spanwise  location  when  the  local  value  of  LESP  equals  the  two-dimensional  value  of 
critical  LESP.  Comparison  of  the  low-order  prediction  of  LEV  initiation  for  a  large  number  of  wing  planforms 
with  CFD  predictions  is  excellent,  confirming  that  the  flow  physics  governing  LEV  initiation  in  finite  wings  is 
the  same  as  that  for  airfoils.  Further,  the  critical  value  of  LESP  obtained  for  two-dimensional  airfoil  flow  for 
one  "high-rate"  motion  can  be  used  for  prediction  of  LEV  initiation  on  any  finite-wing  geometry  undergoing 
any  other  "high-rate"  motion.  Discrepancies  between  the  low-  and  high-order  prediction  of  LEV  initiation 
were  seen  in  very  low-aspect-ratio  wing  (AR  =  2),  for  which  the  rolled-up  shed  vorticity  from  the  wing  tip  has 
a  significant  influence  on  the  wing  flow.  Since  the  UVLM  does  not  model  the  tip-vortex  roll  up,  there  is  a 
noticeable  discrepancy  for  the  AR  =  2  geometry. 

The  second  part  of  the  research  focused  on  modifying  the  UVLM  to  handle  LEV  formation.  In  this  modified 
UVLM,  a  vortex  sheet  from  the  leading  edge  is  modeled,  and  its  geometry  is  calculated  in  each  time  step  by 
assuming  that  the  corner  points  of  the  panels  forming  the  sheet  convect  with  the  local  velocity.  The 
resulting  self-induced  roll  up  of  the  sheet  is  simulated,  with  the  geometry  agreeing  reasonably  well  with 
CFD  results.  While  the  low-order  prediction  of  LEV  formation  shows  promise,  significant  challenges  remain 
in  the  modeling  of  the  sheet  roll  up  and  its  intersection  with  the  wing  geometry.  Discrete-vortex 
amalgamation  techniques,  studied  as  part  of  the  current  effort  for  2D  flows,  shows  some  promise  for 
reducing  the  number  of  vortex  lattices  in  LEV  sheets  from  finite  wings.  It  is  likely  that  some  adaption  of 
vortex  amalgamation  applied  to  vortex  sheets  from  leading  edges  will  allow  for  robust  modeling  for  LEV 
shedding  on  finite  wings. 

Suggested  future  work  includes  (i)  augmentation  of  the  current  UVLM  formulation  to  include  better 
modeling  of  tip-vortex  shedding,  (ii)  development  of  a  hybrid  vortex-sheet /vortex-filament  model  for  LEV 
shedding  from  f i '[i)i Jfnbutio^  piWn'c^ lease01 1  e  LEV  formation,  growth,  and 


detachment,  allowing  for  modeling  of  intermittent  LEV  shedding. 
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